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We present a first simulation of the post-merger evolution of a black hole-neutron star binary in full general 
relativity using an energy-integrated general relativistic truncated moment formalism for neutrino transport. We 
describe our implementation of the moment formalism and important tests of our code, before studying the 
formation phase of an accretion disk after a black hole-neutron star merger. We use as initial data an existing 
general relativistic simulation of the merger of a neutron star of mass 1.4M© with a black hole of mass 7M© and 
dimensionless spin xbh = 0.8. Comparing with a simpler leakage scheme for the treatment of the neutrinos, we 
find noticeable differences in the neutron to proton ratio in and around the disk, and in the neutrino luminosity. 

We find that the electron neutrino luminosity is much lower in the transport simulations, and that both the disk 
and the disk outflows are less neutron-rich. The spatial distribution of the neutrinos is significantly affected by 
relativistic effects, due to large velocities and curvature in the regions of strongest emission. Over the short 
timescale evolved, we do not observe purely neutrino-driven outflows. However, a small amount of material 
(3 x 10 _4 M 0 ) is ejected in the polar region during the circularization of the disk. Most of that material is 
ejected early in the formation of the disk, and is fairly neutron rich (electron fraction Y e ~ 0.15 — 0.25). 

Through r-process nucleosynthesis, that material should produce high-opacity lanthanides in the polar region, 
and could thus affect the lightcurve of radioactively powered electromagnetic transients. We also show that by 
the end of the simulation, while the bulk of the disk remains neutron-rich (Y e ~ 0.15 — 0.2 and decreasing), 
its outer layers have a higher electron fraction: 10% of the remaining mass has Y e > 0.3. As that material 
would be the first to be unbound by disk outflows on longer timescales, and as composition evolution is slower 
at later times, the changes in Y e experienced during the formation phase of the disk could have an impact on 
nucleosynthesis outputs from neutrino-driven and viscously-driven outflows. Finally, we find that the effective 
viscosity due to momentum transport by neutrinos is unlikely to have a strong effect on the growth of the 
magnetorotational instability in the post-merger accretion disk. 

PACS numbers: 04.25.dg, 04.40.Dg, 26.30.Hj, 98.70.-f 


I. INTRODUCTION 

The likely detection of gravitational waves by the advanced 
LIGO/VIRGO detectors (HE 1 in the coming years will open 
up an entirely new way to observe the universe, complement¬ 
ing existing electromagnetic and neutrino observations. Merg¬ 
ers of black holes and neutron stars are expected to be among 
the first and most common sources of gravitational waves to 
be observed 0. Beyond the excitement associated with the 
first gravitational wave detection, binary mergers will provide 
us with a wealth of information which could constrain the 
formation and evolution of massive binaries ala, the out¬ 
comes of core collapse supernovae, or the properties of the 
cold, dense neutron-rich matter in the core of neutron stars El¬ 
ia. 

In the presence of at least one neutron star, the gravitational 
wave signal may be accompanied by electromagnetic and neu¬ 
trino emissions. The joint detection of a system by a grav¬ 
itational wave observatory and an electromagnetic telescope 
would provide additional information about the location, en¬ 
vironment, and parameters of the binary Col. It could also 


help us understand some important astrophysical processes. 
Neutron star mergers are among the most likely progenitors 
of short-hard gamma-ray bursts. They may also significantly 
contribute to the production of heavy elements in the universe, 
through r-process nucleosynthesis in the neutron-rich material 
unbound during some mergers ehhi. This nucleosynthesis 
could be observed through radioactively powered transients, 
observable in the optical and/or infrared bands days after the 
merger (“kilonovae”) [14 , 15 j. 

To understand which signals can be emitted by a merger, 
and how they depend on the initial parameters of the binary, 
we need numerical simulations in full general relativity: the 
strong non-linearities in the evolution of the spacetime sur¬ 
rounding the binary make any approximate treatment of grav¬ 
ity unreliable. Most general relativistic simulations focus on 
a very short period around the time of merger 100 ms for 
simulations involving neutron stars), when general relativistic 
effects are important. Over that time period and beyond grav¬ 
ity, the most important physical effects to take into account are 
the equation of state of neutron star matter, the effects of mag¬ 
netic fields, and the cooling and composition evolution due to 
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neutrino-matter interactions. Although numerical simulations 
of these effects have significantly improved over the last few 
years, much work is needed to simulate them well-enough to 
reliably predict the post-merger signals produced by compact 
binary mergers. The main object of this work is the develop¬ 
ment and testing of an improved method for the treatment of 
neutrinos in the SpEC code COD, based on the moment for¬ 
malism of Thorne 03 . and its application to the post-merger 
evolution of a black hole-neutron star binary. 

In terms of the neutron star equations of state, most general 
relativistic simulations use either a simple Gamma-law equa¬ 
tion of state, or parametrized piecewise-polytropic equations 
of state with a thermal Gamma-law. When modeling the grav¬ 
itational wave signal up to the point of merger, this should be 
sufficient: the gravitational wave signal mostly probes a sin¬ 
gle parameter of the equation of state, its tidal deformability, 
while the detailed structure of the neutron star appears unim¬ 
portant □ El. But the equation of state plays a crucial role 
when simulating the merger and post-merger evolution of the 
binary ini ei. Only a few simulations have used cold cir 
ED or hot 03, [224751 tabulated, nuclear-theory based equa¬ 
tions of state, and of those only the equation of state used 
in na appears consistent with the most recent models for 
dense neutron-rich matter [ 26, 271- I n this work, we use a 
hot, composition dependent, nuclear-theory based equation of 
state by Lattimer & Swesty [28 ], which leads to an acceptable 
mass-radius relationship for neutron stars and allows us to in¬ 
clude neutrino-matter interactions, but is not fully consistent 
with the latest constraints from nuclear experiments. 

The simulation of magnetic fields has seen some impres¬ 
sive recent improvements. General relativistic simulations of 
binary neutron star mergers assuming ideal magnetohydrody¬ 
namics have been performed with a number of codes (see 
e.g. 1291 for a recent review of the field), but resolving the 
growth of magnetic instabilities at an acceptable computa¬ 
tional cost remains an open problem lf30l . Simulations us¬ 
ing force-free or resistive magnetohydrodynamics, needed to 
properly simulate magnetically dominated regions, have also 
been performed [ 3j}j33]. So far they have, however, mostly 
studied the pre-merger evolution of the system. Here, we 
do not take the magnetic fields into account at all, and focus 
solely on neutrino effects. 

The third component needed to simulate compact binaries 
around the time of merger, neutrinos, is also the most com¬ 
putationally expensive to treat accurately. In theory, for each 
species of neutrinos, one should evolve the distribution func¬ 
tion of neutrinos in both physical and momentum space. The 
high dimensionality of the problem places it out of reach of 
current numerical codes and computers. Accordingly, various 
approximations have been developed - most of them in New¬ 
tonian simulations and/or with the aim of studying core col¬ 
lapse supernovae. In general relativistic simulations of binary 
mergers, a few simulations have used a simple cooling pre¬ 
scription (“leakage”) to model the first-order impact of neutri¬ 
nos on the cooling of the disk and the composition evolution 
of high-density regions 1 2314251 . The only published result us¬ 
ing a more advanced method ED studied a binary neutron star 
merger using an energy-integrated (“gray”) version of the mo¬ 


ment formalism for radiation transport ifTTl l34l l35l . The first 
two moments of the neutrino distribution function (i.e. the 
energy density and flux of neutrinos) were evolved, and the 
analytical Ml closure was used to compute the third moment 
(pressure tensor). In this work, we discuss the implementation 
of a similar scheme within the SpEC code, and its applica¬ 
tion to the post-merger evolution of a black hole-neutron star 
merger. As in CD , we use a gray Ml scheme. We should note 
that there are significant differences between our implementa¬ 
tion of the Ml formalism and the one used i n fl3l . in partic¬ 
ular in the treatment of optically thick region^] which would 
make a comparison between the results of the two codes an 
interesting test of the accuracy of the Ml formalism in the 
core of the disk or in the presence of a hypermassive neutron 
star. However, we will not discuss this here, as we are consid¬ 
ering a completely different physical setup for the evolution, 
and the details of the Ml implementation used in [ED have 
not yet been published. 

The focus of this paper is to discuss three different ques¬ 
tions. First, we describe the implementation and testing of 
the Ml formalism in SpEC. An overview of the methods is 
offered in Sec. [ll| while the interested reader will find the de¬ 
tails of the algorithm and important tests of the code in the 
Appendices. Second, we simulate the evolution of a black 
hole-neutron star binary in the critical phase in which the ac¬ 
cretion disk is formed, from the moment at which neutrino 
emission increases due to the heating of the forming disk, to 
the time at which the disk reaches a quasi-equilibrium state 
(about 20 ms after merger). We pay particular attention to the 
temperature and composition evolution, the ejection of matter 
along the spin axis of the black hole, the geometry of the neu¬ 
trino radiation, and the potential astrophysical consequences 
of our results. This is the main focus of Secs. |ry]Tv| And third, 
we study the impact of using the Ml formalism instead of the 
simpler leakage scheme, and the impact of the various approx¬ 
imations which have to be made when using a gray scheme. 
This is the focus of Sec. lIV Gl 


II. MOMENT FORMALISM 

We will start here with an overview of the Ml formalism, 
and of its numerical implementation in SpEC. A more detailed 
description of many aspects of the algorithm, as well as tests 
of our implementation, are available in the appendices. 


A. Evolution equations 

For each neutrino species z^, we can describe the neutri¬ 
nos by their distribution function /( Z y)(^ a ,p a ), where x a = 
(t, x l ) gives the time and the position of the neutrinos, and p a 
is the 4-momentum of the neutrinos. The distribution function 
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/(V) evolves according to the Boltzmann equation 
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where the T^ are the Christoffel symbols and the right-hand 
side includes all collisional processes (emissions, absorptions, 
scatterings). In general, this is a 7-dimensional problem which 
is extremely expensive to solve numerically. Approximations 
to the Boltzmann equation have thus been developed for nu¬ 
merical applications. In this work, we consider the moment 
formalism developed by Thorne fl7lL in which only the low¬ 
est moments of the distribution function in momentum space 
are evolved. Our code is largely inspired by the implemen¬ 
tation of Thorne’s formalism into general relativistic hydro¬ 
dynamics simulations proposed by Shibata et al. m and 
Cardall et al. ESI. We limit ourselves to the use of this for¬ 
malism in the “gray” approximation, that is we only consider 
energy-integrated moments. Although the moment formalism 
can in theory be used with a discretization in neutrino en¬ 
ergies, this makes the simulations significantly more expen¬ 
sive and involves additional technical difficulties in the treat¬ 
ment of the gravitational and velocity redshifts. We will also 
only consider three independent neutrino species: the elec¬ 
tron neutrinos z/ e , the electron antineutrinos z/ e , and the heavy- 
lepton neutrinos v x . The latter is the combination of 4 species 
(z/^, v T , z/ r ). This merging is justified because the temper¬ 
atures and neutrino energies reached in our merger calcula¬ 
tions are low enough to suppress the formation of the corre¬ 
sponding heavy leptons whose presence would require includ¬ 
ing the charged current neutrino interactions that differentiate 
between these individual species. 

In the gray approximation, and considering only the first 
two moments of the distribution function, we evolve for each 
species projections of the stress-energy tensor of the neutrino 
radiation T^ d . One possible decomposition of T^ d is l34l 

T rad = + H»U V + H V U» + S pV , (2) 


with H p Up = = 0 and u p the 4-velocity of the fluid. 

The energy J, flux H p and stress tensor S pL/ of the neutrino 
radiation as observed by an observer comoving with the fluid 
are related to the neutrino distribution function by 


j = 

pOO 

L dvv3 j 

[ <Klf {u) (x a : 

, v , 0) , 

(3) 

H>' = 

pOO 

L ivv3 i 

[ (x a . 


(4) 

S^ = 

poo 

L 

[ dflf( v )(x a . 


(5) 


where v is the neutrino energy in the fluid frame, f dfl de¬ 
notes integrals over solid angle on a unit sphere in momentum 
space, and 


p a = v(u a + l a ) , 


( 6 ) 


with l a u a = 0 and l a l a = 1. We also consider the decom¬ 
position of T^ d in terms of the energy, flux and stress tensor 
observed by an inertial observer, 

T rad = En ^ nU + Fl±nV + + FlXV , (7) 


with = P pv np = F 1 = P tv = 0, and n a the unit 

normal to a t = constant slice. 

We use the 3+1 decomposition of the metric, 

ds 2 = g a pdx a dx& ( 8 ) 

= —a 2 dt 2 + 7 ij(dx l + /3 l )(dx J + /J* 7 ') , (9) 

where a is the lapse, fd l the shift, and 7 ^ the 3-metric on a 
slice of constant coordinate t. The extension of 7 ^ to the full 
4-dimensional space is the projection operator 

7 = Qa(3 + n a np . (10) 

We similarly define a projection operator onto the reference 
frame of an observer comoving with the fluid, 

hat(3 = 9a/3 "T U a Uf3 • (H) 

We can then write equations relating the fluid-frame variables 
to the inertial frame variables lf35l : 

E = W 2 J + 2WvpH p + v^VyS^ , 

= W 2 v /l J + W(g^ - 

+Wv^v u H u + (g^y - n ll v v )VpS up , 

F ^v = W 2 VpV u J + W(gp P - npV p )v„H p 

-\~W (g pv pF[ p 

TipV p)(s)vk TivVfd)S P , 

using the decomposition of the 4-velocity 

u p — W(n p + v p ) , 

with = 0 and W = ^/l + 7 ^muj. 

Evolution equations for E = E and F = ^ 7 F 1 can 
then be written in conservative form: 

dtE + djiaF* -0 j E) (16) 

= a(P^Kij - F j dj In a - S a n a ) , 
dtFi + djiaP? -ftFi) (17) 

= (-Edict + F k di/3 k + ^P jk 8 iljk + aS a 7ia ) , 

where 7 is the determinant of 7 ^, P{j = ypyPij , and S a in¬ 
cludes all collisional source terms. 

To close this system of equations, we need two addi¬ 
tional ingredients: a prescription for the computation of 
P^(E,Fi) (‘closure relation’, which we choose following 
Minerbo 1978 l36l k and the collisional source terms S a . In 
the Ml formalism, the neutrino pressure tensor P l i is recov¬ 
ered as an interpolation between its known limits for an op¬ 
tically thick medium and an optically thin medium with a 
unique direction of propagation for the neutrinos. We pro¬ 
vide details on its computation in Appendix [A| For the source 
terms S a , we will consider that the fluid has an emissivity rj 
due to the charged-current reactions 

p + e~ -+ n + v e , 
n + e + -+ p + v e , 


( 12 ) 

(13) 

(14) 

(15) 


(18) 

(19) 
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as well as electron-positron pair annihilation 

e + + e~ ViVi , (20) 

plasmon decay 

7 -» Wi , ( 21 ) 

and nucleon-nucleon Bremsstrahlung 

+ + 27 + 27 . ( 22 ) 

The inverse reactions are responsible for an absorption opacity 
K a . We also consider a scattering opacity k s due to elastic 
scattering of neutrinos on nucleons and heavy nuclei. The 
source terms are then 

S a = ^y(r]u a - K a Ju a - (K a +K s )H a ) . (23) 

Details of the choices made for the computation of the gray rj, 
K a and k s are provided in Appendix [B] 

B. Numerical scheme 

We add the evolution of neutrinos with the moment scheme 
to the SpEC code lfl 6 l . which already includes a general rela¬ 
tivistic hydrodynamics module l37l . The latest methods used 
for evolving in SpEC the coupled system formed by Einstein’s 
equation and the general relativistic equations of hydrody¬ 
namics are described in 1381 , Appendix A. 

In the Ml formalism, we evolve the variables E and Fi 
according to Eqs. ( [T6|[T7] ), and couple the neutrinos to the fluid 
evolution. The coupling takes the form 

d t r = ... + aS a n a , (24) 

dtSi = ... , (25) 

d t {p*Y e ) = ... - sign(Pi)a^ ~ - Kc — , (26) 

\ e W 

where p*, p*Y e ,r and Si are the conservative hydrodynamics 
variables which are evolved in SpEC, 


where 

roo k 

F M = ) dx (32) 

Jo 1 n- exp [x r/jy) 

is the Fermi integral, rj u = p u /T is the degeneracy parameter, 
and fi v is the chemical potential of neutrinos in equilibrium 
with the fluid, (e 2 leak ) is the global estimate of the average 
square energy of neutrinos obtained from the simpler leak¬ 
age scheme l24l . and (e 2 fi uid ) would be the average square 
energy of the neutrinos if they obeyed a Fermi-Dirac distribu¬ 
tion at the local temperature T of the fluid and with neutrino 
chemical potential fi v (see Appendix [b}. This form for (e u ) 
is exact when n a oc e 2 , the neutrinos are in equilibrium with 
the fluid, and (el leak ) < (el Jluid ) - conditions which, for v e 
and z/ e , are approximately satisfied in optically thick regions. 
In particular, the last condition is satisfied in post-merger ac¬ 
cretion disks because {e 2 leak ) is set by the temperature on 
the neutrinosphere, which is generally colder than points of 
higher optical depth in the disk. When one of these conditions 
is not satisfied, on the other hand , it is one of the approxi¬ 
mations which we have to make to accommodate the use of a 
gray scheme, one that is well motivated by the v 2 dependence 
of the dominant absorption processes, and the relative homo¬ 
geneity of the fluid temperature around the neutrinospheres. 

For the applications that we are considering here, the back- 
reaction of the neutrinos onto the fluid evolution is weak (ex¬ 
cept for transients when we turn on neutrino emission). Ac¬ 
cordingly, we separate the hydrodynamics and neutrino evo¬ 
lution. Our evolution scheme thus proceeds as follow: 

• Evolve the Einstein equations and the general relativis¬ 
tic hydrodynamics equations, without taking neutrinos 
into account, over a time step A tn chosen as in [38), 
Appendix A.3. 

• Evolve the neutrino radiation, potentially taking mul¬ 
tiple time steps, so that the neutrinos are also evolved 
by A tn- The time step used to evolve the neutrinos 
is chosen as described in Appendix [D] Each time step 
proceeds as follow: 


P * = PoW^f , (27) 

r = p*{hW - 1) - v^P , (28) 

Si = p* hui , (29) 

po is the baryon density of the fluid, P its pressure, Y e its elec¬ 
tron fraction, and h its specific enthalpy. (e u ) is the weighted 
average energy of neutrinos, which should be 

/ V _ fo°° - Ka(ev)J(e„)\ de u 

[eu) roo , > IJUj 

Jo e u atl/ 

and sign(^) is 1 for is e , —1 for 9 e , and 0 for heavy-lepton 
neutrinos. Lacking the knowledge of the neutrino spectrum, 
we use for (e u ) the approximation 


(tv) 


F 4 M 


T max 



(31) 


- Reconstruct the fields E,Fi/E at cell faces, using 
the minmod reconstruction method. 


- Compute the closure relation at cell faces to get 
the fluxes ( aF l — /PP) and (aPj — 


-Use those fluxes to compute the divergence terms 
in Eqs. using the shock-capturing meth¬ 

ods and corrections in high optical depth regions 
described in Appendix [C] 


- Compute the closure relation at cell centers. 


- Compute the gravitational source terms on the 
right-hand side of Eqs. ( pPp7| ) [everything but the 
terms proportional to S a ] from E and F, at the be¬ 
ginning of the neutrino step. 


- Solve Eqs. ( I6p7] ) by treating implicitly the colli- 
sional source terms proportional to S a , following 
the method described in Appendix [D] 
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- Compute the coupling to the hydrodynamics 
variables, and update (r, Si,p*Y e ) according to 
Eqs. ([24]|26]). 

- If we have evolved the neutrinos by A/#, go back 
to the GR-Hydro evolution. Otherwise, take the 
next neutrino time step. 

A more complete description of the different steps of this 
algorithm is provided in Appendices |A||D| 

III. INITIAL CONDITIONS AND NUMERICAL SETUP 

As a first astrophysical application of our code, we consider 
the disk formation phase of a black hole-neutron star merger 
from Foucart et al. 2014 (241. This phase is particularly in¬ 
teresting to study with a general relativistic code and neutrino 
transport because general relativistic effects and the evolution 
of the metric remain important at this point. Furthermore due 
to the high temperatures experienced by the fluid during disk 
formation, the neutrino luminosity is higher and the composi¬ 
tion evolution is faster than at any other time. Finally the fluid 
is initially very close to the black hole, where relativistic ef¬ 
fects cannot be neglected, and it is far from equilibrium. This 
phase of the evolution would thus not be properly captured by 
simulations which start from an equilibrium torus configura¬ 
tion, or which model general relativistic effects through the 
use of pseudo-Newtonian potentials. 

In the specific merger that we are considering, the masses 
of the compact objects before merger are M^ s = 1.4M e and 
Mg H = 7M 0 . The initial dimensionless spin of the black 
hole is Xbu = 0 - 8 , and it is aligned with the orbital angular 
momentum of the binary. The neutron star is initially non¬ 
spinning. We showed in f24l that the disruption of the neutron 
star results in the ejection of about 0.06M o of material, and 
the formation of an accretion disk of mass Mdi s k ~ 0.1 M 0 . 
The final properties of the black hole are Mg H = 8 M 0 and 
Xbh = 0-87. We use the equation of state of Lattimer & 
Swesty [28 ] with nuclear incompressibility parameter Kq = 
220 MeV and symmetry energy S„ = 29.3 MeV, using a table 
available on http://www.stellarcollapse.org and described in 
O’Connor & Ott 2010 |[39l . For a neutron star of mass M N s = 
1.4M 0 , this results in a neutron star radius 77 ns = 12.7 km. 

The post-merger configuration obtained as a result of the 
merger is expected to be fairly typical for black hole-neutron 
star mergers in which the neutron star is disrupted by the 
black hole. In ED, we studied the range of initial black 
hole masses Mg H = (7 — 10) M 0 and neutron stars masses 
M^g = (1.2 — 1.4) M 0 currently deemed most likely from 
the observation of galactic stellar mass black holes f4(£ 411 
and of neutron stars in compact binary systems |j42l-44jl. The 
distribution of neutron star and black hole masses may be dif¬ 
ferent for black hole-neutron star binaries, but no such binary 
has been observed so far. For those parameters, a moderate to 
high initial black hole spin, %bh — (0-5 — 0.9), is a require¬ 
ment for the neutron star to be disrupted by the black hole 
and thus allow the formation of an accretion disk EH. Once 
that condition is satisfied, however, we showed that the local 


properties of the disk are largely independent of the binary pa¬ 
rameters. 10 ms after merger, we observe T ~ 5 — 10 MeV, 
Y e rsj 0.15 and a typical density p™* ~ 10 10-11 g/cm 3 . 
Global properties, such as the neutrino luminosity, mostly 
scale with the mass of the disk. The accretion disks have 
masses Mdisk ^ (0.05 — 0.15) M 0 , optical depths of a 
few ( r Ve ~ 5, r^ e ~ r Vx ~ 1), and neutrino luminosities 
Ly e ~ Ly e ^ 10 53 erg/s and L„ x ~ 10 52 erg/s for all 4 
species combined. The only exceptions appear to be the ac¬ 
cretion disks formed in mergers with rapidly rotating, rela¬ 
tively low mass black holes (see also ED HD), which can be 
much more massive [Mdisk ~ (0.2—0.5) M 0 ], more optically 
thick, and more luminous (particularly in heavy lepton neutri¬ 
nos v x ). They may also develop additional hydrodynamics 
instabilities causing long-term heating in the disk lf23l . 

In this paper, we want to assess the effects of a better treat¬ 
ment of the neutrinos on the evolution of the disk and of its 
immediate neighborhood. To do so, we compare simulations 
using the same leakage scheme as in (23j|24l with simulations 
using the Ml formalism presented here. The leakage scheme 
should give us a rough estimate of the cooling of the disk due 
to the neutrinos and of the evolution of the composition of the 
high density regions of the disk, but does not include heat¬ 
ing or composition evolution due to neutrino absorption in 
the low-density regions. The Ml formalism should give us 
a better estimate of both the total luminosity and the effects of 
neutrino-matter interactions. 

We evolve the post-merger accretion disk with the follow¬ 
ing treatments of the neutrinos: 

• A “Leakage” simulation using the same algorithm as 
in ESJ EO . 

• A “Leakage” simulation in which the neutrino opacities 
have been corrected as in our Ml code (to guarantee 
that the energy density of neutrinos in the optically thick 
regions is correct, see Appendix [B]). 

• A “Ml” simulation using the standard methods de¬ 
scribed in this paper. 

• A “Ml” simulation with a simpler treatment of the 
emissivity in optically thin regions (i.e. without the cor¬ 
rection described in Eq. |B 6 | ), to get a first estimate of the 
errors due to the use of a gray scheme: we find in our 
test of the neutrino-fluid interactions in a post-bounce 
supernova profile (see Appendix |E 6 | ) that variations in 
the results with and without that correction are compa¬ 
rable to the difference between the gray results and the 
results of an energy-dependent code. 

• Three “Ml” simulations with a numerical grid covering 
a larger volume and improved boundary conditions (see 
below), with and without the correction to the emissiv- 
ities described in Eq. |B 6 [ and with and without the cor¬ 
rection to the opacities described in Eq. |B 8 | 

As initial configuration, we consider a snapshot of simula¬ 
tion M14-7-S8 of (241 at to = t merge + 6.1 ms, where t merge 
is the time at which half of the neutron star material has been 
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FIG. 1. Density and velocity field in the equatorial plane of the black 
hole at the initial time to ~ /merge + 6.1ms. For reference, the 
velocity at the peak of the density distribution is about 0.6c. 


accreted by the black hole. This is around the time when an 
accretion disk is forming and neutrino effects begin to signif¬ 
icantly affect its evolution. The initial configuration is shown 
in Fig. [I] From the disruption of the neutron star to to, the 
hydrodynamics equations were evolved using two levels of 
refinement, each with 100 2 x 50 grid points and with grid 
spacings Ax « (1.5 km, 3 km) in the equatorial plane and 
Az ~ (0.5 km, 1 km) in the vertical direction. For compari¬ 
son, the horizon of the black hole has a radius vbu ~ 25 km in 
the coordinates of our simulation, the peak of the matter den¬ 
sity in the accretion disk is initially at rui s k ^ 50 km, and the 
scale height of the disk is H ~ (0.2 — 0.3)r. We chose the grid 
structure in order for the finest grid level to cover the forming 
accretion disk. We also used the equatorial symmetry of the 
system to only evolve the 2 > 0 region. Note that, in SpEC, 
Einstein’s equations are evolved on a separate pseudo spectral 
grid extending 3700 km away from the center of the black 
hole. This is why we can use such a small finite difference 
grid if we only want to study the evolution of the post-merger 
accretion disk. A smaller grid means that we can maintain 
a reasonable resolution in the disk at a fairly low computa¬ 
tional cos|5 but also that the unbound material, and some of 
the bound material in the tidal tail formed when the neutron 
star is disrupted, is allowed to leave the grid. By to, we still 
have a baryonic mass M\b = 0.12M 0 on the grid, while we 
have allowed O.O6M 0 of unbound material and 0.06M 0 of 
bound material to escape. In (24), we also performed lower 
resolution simulations following the material in the tidal tail 
farther out. From these simulations, we can deduce that the 
material that we allowed to escape would start to fall back 


2 Although such a resolution is only reasonable because we do not include 

the effects of magnetic fields. It would be much too coarse to capture the 
effects of MHD instabilities. 


onto the disk at t ^ /merge + 20 ms. At later times, material 
neglected in our simulation may affect the evolution of the 
disk. 


Our simulations do not include the effects of magnetic 
fields, and the only viscosity is due to the finite resolution 
of our numerical grid. Evolving for much longer would thus 
not be very informative. The growth of magnetic fields due 
to the magnetorotational instability is expected to act on a 
timescale comparable to the orbital timescale of the disk, 
i.e. a few milliseconds, and angular momentum transport 
due to magnetic turbulence is expe cted to act on timescales 
/vise 00 (30 — 700) ms (see Sec. IV A). Accordingly, we 
only study the evolution of the post-merger disk from to to 
tf = /merge + 20 ms. We will see that this is long enough for 
the fluid and the neutrinos to reach a quasi-equilibrium state. 
This also allows us to study the main differences between sim¬ 
ulations using the leakage scheme and simulations using the 
Ml formalism. Finally, we can get a first estimate of the out¬ 
flows emitted in the region above the disk and of the effects of 
interactions between the disk and the tidal tail. 


For the evolution between to and //, we use a different 
finite difference grid. Indeed, although the system still re¬ 
spects the equatorial symmetry after merger, small perturba¬ 
tions which are not equatorially symmetric might grow in the 
disk due to hydrodynamical instabilities. To make sure that 
our grid structure does not artificially suppress such pertur¬ 
bations, we remove the assumption of equatorial symmetry. 
We also want to capture the radial growth of the disk, poten¬ 
tial matter outflows, and the evolution of the bound material 
in the tidal tail remaining on the grid (which is still expand¬ 
ing at / = to). For the post-merger evolutions, we thus use 
a third level of refinement (with Ax ~ 6 km in the equato¬ 
rial plane and Az ~ 2 km in the vertical direction), keeping 
as before Ax « (1.5 km, 3 km) in the equatorial plane and 
Az ~ (0.5 km, 1km) in the vertical direction for the finer 
levels. Each refinement level has 100 2 3 grid points. In three of 
the Ml simulations, a fourth level is added with Ax = 12 km 
and Az = 4km. This was done because, with only 3 levels 
of refinement, we cannot follow the evolution of low-density 
material far enough along the direction of the black hole spin 
axis to extract information about potential disk winds. We 
note that, due to our choice of map between the grid and the 
laboratory frame, the resolution in the grid is not uniform. We 
quote here the resolution close to the black hole, but the grid 
spacing far from the black hole is actually ^30% coarser. 

We also performed simulations with different numerical 
resolution, to assess our numerical errors. A detailed discus¬ 
sion of these errors is provided in Appendix |F| In summary, 
we find that the errors due to finite resolution are at most com¬ 
parable to the errors due to the use of a gray scheme, and much 
smaller than the errors in the leakage scheme. 

As opposed to previous SpEC simulations, we want here 
to study the behavior of low-density regions. Accordingly, for 
these last Ml simulations we also made a few modifications to 
the handling of low-density matter. In general relativistic hy¬ 
drodynamics simulations, the region surrounding the neutron 
star or accretion disk is generally filled with lower-density ma¬ 
terial in order to avoid evolving towards either negative den- 
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sities or unphysical values of the evolved variables. Addition¬ 
ally, corrections are applied to all regions of baryon density 
po lower than an “atmosphere” threshold Po tm . In our case, 
those corrections set the temperature to T = 0.5 MeV and 
the spatial components of the 4- velocity to Ui = 0, although 
different prescriptions also work in practice. This generally 
results in large portions of the grid being filled with mate¬ 
rial at po ~ pg tm . In previous SpEC simulations with the 
LS220 equation of state 1231 l24l . that threshold was set at 
Po tm = lO -5 ^™^, with the maximum value of the 

baryon density on the grid at the current time. In accretion 
disks, this leads to po tm ~ 10 6_7 g/cm 3 * . Due to larger nu¬ 
merical errors close to the black hole, Po tm i s a l so higher in a 
region about one apparent horizon radius away from the black 
hole horizon, by as much as a factor of 100 on the horizon 
itself. In this work, we are interested in disk outflows with 
density po ^ 10 7_8 g/cm 3 in the regions covered by our 
grid. Such outflows cannot be launched with that high an at¬ 
mosphere threshold. This is in part because of the pressure 
exerted by the atmosphere material, and in part because any 
material decompressing to Po tm is immediately slowed down 
to Ui = 0, and falls back onto the disk. We thus modify our 
choice of atmosphere threshold to 


_atm _ .floor 

P 0 — P 0 


10 3 exp 



(33) 


with tah the grid-coodinate radius of the black hole apparent 
horizon, which is by construction a sphere in the coordinates 
of our numerical grid, r is the distance to the center of the 
black hole in the same coordinates, and Po° or ^ 10 5 g/cm 3 . 
With this choice, the atmosphere threshold remains as before 
close to the black hole, but now drops to po tm ~ Po°° r at 
larger distances, which is sufficient for our current purpose. A 
lower threshold would be required if we wanted to follow disk 
outflows farther from the black hole. 

Finally, in previous simulations the table containing 
the equation of state information only covered the range 
10 8 g/cm 3 < po < 10 16 g/cm 3 , and was artificially extended 
to lower densities using a simple ideal gas equation of state. 
In this last simulation, we instead use a table going down to 
10 5 g/cm 3 . The new table uses 345 logarithmically spaced 
points to cover that range of densities, instead of 250 for the 
smaller table. Both tables are also discretized in temperature 
(0.01 MeV < T < 251 MeV with 136 logarithmically spaced 
points) and composition (0.035 < Y e < 0.53 with 50 lin¬ 
early spaced points). Below 10 8 g/cm 3 , the compositional 
information of the equation of state is approximated by the 
nuclear statistical equilibrium of matter at 10 8 g/cm 3 , for the 
same temperature and electron fraction (see O’Connor & Ott 
2010 l39l for details). 


IV. DISK EVOLUTION 

In the following, we describe the physical results of the 
simulations. We mainly focus on the Ml simulation using 
our standard algorithm, including the improved treatment of 
the low-density regions for the hydrodynamical variables, as 
it offers the most accurate predictions. The effects of different 
neutrino treatments are discussed in Sec. II V Gl 


A. Global properties of the disk 

At the beginning of the simulation (to = t merge + 6.1 ms), 
we have about O.O8M 0 of material in a forming accretion disk 
extending ~ 60 km away from the black hole. As shown in 
Fig. [T] this material is still far from being circularized or ax- 
isymmetric. Another 0.04M 0 of material is on the grid in 
an extended tidal tail with a relatively flat density profile. As 
discussed in the previous section, an additional 0.06 M 0 of 
unbound material and 0.06 M 0 of bound material with fall¬ 
back time longer than ^ 20 ms were allowed to escape the 
grid before to- The initial configuration has a sharp density 
peak around r = 45 km. Most of the tail material is cold 
(T < 1.5 MeV), while the disk material has a broad tem¬ 
perature distribution, with most of the matter at temperatures 
2 MeV < T < 12 MeV. The composition of the disk and tail 
is sharply peaked at 0.05 < Y e < 0.07, and the small amount 
of material at Y e > 0.1 is in the hot regions close to the black 
hole and rapidly accreted. 

Although neutrinos have an important effect on the evolu¬ 
tion of the system, to first order, purely hydrodynamical ef¬ 
fects dictate the evolution. During the time period considered 
in this simulation, 5 ms — 20 ms after merger, the circular¬ 
ization of the disk material is the most important effect. At 
the initial time, most of the disk material is at or close to pe- 
riastron. Shock heating and the contraction of the disrupted 
material cause the fluid to heat to its maximum average tem¬ 
perature, (T) = 6.4 MeV at t = to + 1.5 ms[^] Afterwards, 
the forming disk goes through a damped cycle of expansions 
and contractions, with a period of about 6 ms. During each ex¬ 
pansion period, the temperature of the fluid and the accretion 
rate decrease. During contractions, shock heating causes the 
temperature to rise, and the accretion rate increases. Neutrino 
emissions and absorptions, although not critical to the dynam¬ 
ics, contribute to a smoothing of the temperature distribution 
and determine the composition of the fluid. Their total lumi¬ 
nosity mostly follows the oscillations of the fluid temperature. 
Energy lost to neutrino emissions and to the accretion of hot 
material onto the black hole also causes a slower global cool¬ 
ing of the fluid, by about 1 MeV over the 14 ms of evolution. 

We look at the evolution of the disk through snapshots of 
the simulation taken at to, t\ = to + 5 ms, t<i = to + 10 ms, 
and tf = to + 14 ms. We plot the fraction of the local mass 


3 Here and in the rest of the text, the average fluid properties refer to density- 

weighted averages 
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FIG. 2. Radial distribution of the material outside of the black hole. 
We plot the fraction of the total mass within bins of width A r ~ 
3 km at 4 representative times of the simulation. Most of the matter 
is initially near the black hole (r < 60 km) and close to periastron. 
The disk then expands, and later circularizes. 



T (MeV) 

FIG. 3. Fraction of the total mass within temperature bins of width 
AT = 0.5 MeV, plotted at 4 representative times of the simulation. 
Within 10 ms, the temperature of the disk becomes very homoge¬ 
neous. 


observed at a given radius in Fig. [2| at a given temperature 
in Fig. [3j and at a given electron fraction in Fig. [4] All plots 
are normalized to the total mass on the grid. Due to accretion 
onto the black hole, the total mass decreases from 0.12 M 0 at 
to to 0.097 M e at A, 0.075 M e at t 2 , and 0.069 M 0 at tf. By 
the end of the simulation, nearly all of the remaining material 
is in a circularized accretion disk. We stop the simulations at 
tf = tmerge + 20 ms, when the material which was allowed to 
leave the numerical grid should begin to affect the evolution 
of the disk. 

The density evolution (Fig. [2]) is mostly a consequence of 
the accretion of material at the inner edge of the disk and the 
circularization of disk and tail material. At first, the large ec¬ 
centricity of most of the disrupted material causes the mat¬ 
ter distribution to broaden. This is most visible in the differ¬ 
ence between to and A, which correspond to times close to 
the moment of maximum contraction and maximum expan¬ 
sion of the forming disk. Afterwards, a more slowly evolving 
disk forms with a density peak r ~ 60 km — 70 km, and some 
milder oscillations. Turbulent angular momentum transport 
due to the magnetorotational instability, which is not modeled 
here, should eventually take over and cause a viscous spread¬ 
ing of the disk. This should happen over the viscous timescale 
Adsc ^ (a£ 2 fi) _1 , where a is the standard viscosity param¬ 
eter, S = H/r and H is the scale height of the disk. For the 
expected a ~ 0.01 — 0.1 and the simulation 5 ~ 0.2 — 0.3, we 
get t v i SC ~ (30 — 700) ms > tf — to = 14 ms. We thus expect 
angular momentum transport due to magnetically-driven tur¬ 
bulence to be fairly unimportant on the timescales considered 
here. 

The temperature evolution (Fig. [3]) is largely a result of mix¬ 
ing in the fluid, shock heating during the circularization of the 
disk, and energy transport by neutrinos. Over the first 10 ms 
of evolution, the temperature goes from a wide, nearly flat dis¬ 
tribution with 2 MeV < T < 12 MeV to a more uniform tem¬ 
perature T ~ 4 MeV—6 MeV. Some mixing of fluid elements 


with different temperatures does occur in purely hydrodynam- 
ical simulations, or in simulations using a simple cooling pre¬ 
scription for the neutrinos (leakage runs). But the process is 
significantly more efficient when using a transport scheme, as 
now neutrinos can transport energy from hot parts of the disks 
to cool parts of the disk. The homogenization of the temper¬ 
ature is helped, in both transport and leakage simulations, by 
the steep dependence of the neutrino emissivities on the tem¬ 
perature, which causes the hottest points in the disk to cool 
much more rapidly than their neighbors. But Fig. [3] shows 
that, in addition to the rapid cooling of the hottest points, neu¬ 
trino transport also causes a heating of the cooler points in the 
disk. Variations in the average temperature of the fluid are due 
mostly to the oscillations observed as the disk circularizes. 

As for the density, magnetically-driven turbulence is ex¬ 
pected to affect the evolution of the temperature over long 
time scales. Newtonian simulations have shown that an 
isotropic viscosity can heat the disk to maximum temperatures 
of T max - 3 MeV - 10MeV for a - 0.001 - 0.1 07). For 
the largest viscosities, viscous heating would be relevant to the 
thermal evolution of the disk towards the end of our simula¬ 
tions. For lower viscosities, it would remain largely irrelevant 
until the disk cools down, but magnetically-driven turbulence 
might still hasten the homogenization of the disk temperature. 
As viscous heating is the main process stopping the cooling 
of the disk, it is of course always relevant to the long term 
evolution of the disk. 

Finally, Fig. [4] shows that the disk also reaches its equilib¬ 
rium composition within about 10 ms. By that point the disk 
is in a quasi-equilibrium state in which charged-current reac¬ 
tions no longer cause net changes of the electron fraction of 
the disk. The composition then evolves more slowly as the 
properties of the disk change, adapting nearly instantaneously 
to the disk evolution. At the temperatures and electron frac¬ 
tions observed in most of the disk, matter becomes less neu¬ 
tron rich when the disk heats up, and more neutron rich as it 
cools down. The electron fraction also decreases for denser 
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FIG. 4. Fraction of the total mass within electron fraction bins of 
width A Ye = 0.01, plotted at 4 representative times of the simula¬ 
tion. Within 10 ms, the electron fraction of the disk becomes very 
homogeneous. The tail of higher- Y e material is due to lower density 
regions around the disk. 
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FIG. 5. Density and velocity field in the equatorial plane of the black 
hole at the end of the simulation, 20 ms after merger. For scale, the 
velocities at the peak of the density distribution are v ~ 0.5c. 


material. At times beyond 20 ms past merger, the disk is ex¬ 
pected to cool and the core of the disk will reneutronize, at 
least until viscous heating balances neutrino cooling. How¬ 
ever, this is not necessarily true for the high-F e tail of the 
distribution which corresponds to cooler, lower density points 
whose composition has been significantly affected by neutrino 
irradiation and which evolve fairly slowly at late times. 


B. Final disk configuration 


From the discussion in Sec. IV A we expect that at tf 
the simulated hydrodynamical properties of the accretion 
disk (density, temperature, composition) are probably a good 
representation of the state of the system about 20 ms after 
merger. Magnetically-driven turbulence and magnetically- 
driven winds could affect the results. But so far simulations 
have found that for reasonable initial values of the magnetic 
fields, magnetic outflows only appear more than 100 ms after 
the merger [48], or are not observed at all l49l l50l |^[ As dis¬ 
cussed in the previous section, MRI-driven angular momen¬ 
tum transport is also expected to occur on timescales longer 
than the duration of our simulation. In this section, we thus 
offer a more detailed description of the disk at the end of the 
simulation, assuming that it is a fairly good description of a 
post-merger accretion disk ~ 20 ms after a black hole-neutron 
star merger. 

We first show the density and velocity field in the equatorial 
plane of the black hole (Fig. [5]), and in the vertical plane y = 
Vbyl (Fig. m >. Here, we define the velocity as the “transport 
velocity” v l T = ^u l , which satisfies d t p* + di{p^v l T ) = 0. 


4 The main difference between the two sets of simulations is the inclusion of 

an initial dipole magnetic field outside of the neutron star in 08). 


This is a convenient definition of the velocity for visualization 
purposes, as it directly shows the motion of fluid elements on 
the grid with all general relativistic coordinate effects taken 
into account. From Fig. [5] we infer that the disk is well circu¬ 
larized up to r ~ 100 km. Beyond that radius, we have either 
the last remnants of the tidal tail falling back onto the disk 
(upper left quadrant), or low density material expanding into 
the neighboring medium. Most of the mass is in the inner re¬ 
gion of the disk, at densities po ~ 10 11_12 g/cm 3 . But beyond 
r ~ 70km, the density drops rapidly to po ~ 10 8_9 g/cm 3 at 
the edge of the disk. The velocity in the inner disk is mildly 
relativistic (vt ~ 0.5c at r ~ 60km) which, as we shall see, 
has important effects on the geometry of the neutrino radia¬ 
tion. 

Fig. [6] shows some interesting additional features. The core 
of the disk at the end of our simulation is of moderate geomet¬ 
rical thickness ( H/r ~ 0.2), and still noticeably asymmetric 
(note that the density is plotted on a logarithmic scale). The 
region close to the spin axis of the black hole is free of any 
material, except for the accretion of low density material from 
the disk at height z < 80 km, with po ~ 10 7-8 g/cm 3 . There 
are however resolved outflows coming from the contact re¬ 
gions between the disk and the tail (left side of the plot), at 
densities po ~ 10 8 g/cm 3 . At the grid boundary, this material 
has not reached the escape velocity (it has u t > — 1). But most 
of the outflowing fluid has enough energy to be unbound (i.e. 
it satisfies the condition hu t < —1, thanks to temperatures 
T 2 MeV). We discuss the disk/tail interactions in more de¬ 
tail in Sec. |IV C[ and the unbound material in Sec. |IVD[ For 
now, we simply note that the outflows, although helped by 
neutrino heating and radiation pressure, are mostly a conse¬ 
quence of purely hydrodynamical interactions at the disk/tail 
interface, as can be verified by their existence in the simula¬ 
tions using a leakage scheme. Earlier in the evolution, shocks 
during the circularization power similar outflows, albeit more 
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FIG. 6. Density (in g/cm 3 ) and velocity field in the y — ?/bh plane 
at the end of the simulation, 20 ms after merger. For reference, the 
largest poloidal velocities are about 0.6c. In the core of the disk, the 
velocities are mostly directed out of the plane of the visualization. 
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axisymmetric and closer to the black hole spin axis. 

On the opposite side of the disk, we observe a more homo¬ 
geneous expansion of low-density material. This expansion 
should in practice be stopped once the outflowing disk mate¬ 
rial begins interacting with the tail material neglected in this 
simulation, i.e. the material falling back onto the disk on a 
timescale longer than the duration of the simulation. 

Similar visualizations for the temperature are shown in 
Fig.[7]and Fig. [8] Asymmetric features in the temperature dis¬ 
tribution remain clearly visible despite the smoothing effect 
of the neutrino cooling and heating. Nearly all of the mate¬ 
rial, including the tidal tail, has been heated to T > 1 MeV, 
by a combination of shocks and neutrino absorption. Denser 
material with po > 10 11 g/cm 3 is heated to T > 4 MeV. 
The relatively small variations of the temperature are one of 
the main reasons a gray scheme for the neutrino radiation is 
more reliable for merger simulations than in core collapse su¬ 
pernovae, where large temperature gradients exist and small 
changes in the interactions between the neutrinos and the fluid 
can significantly modify the dynamics of the system. 

Finally, we consider the composition of the fluid, through 
its electron fraction Y e , in Fig. [9] and Fig. [10| Being able to 
reliably predict the electron fraction in low-density regions 
is one of the main advantage of the use of the moment for¬ 
malism over the leakage scheme. We see that the core of the 
disk remains relatively neutron rich (0.15 < Y e < 0.20), but 
material below po ~ 10 10 g/cm 3 is significantly protonized. 
About 30% of the total mass is at Y e > 0.2, and about 10% 
at Y e > 0.3. More importantly, the less neutron-rich material, 
which surrounds the neutron-rich core of the disk, is more 
likely to be unbound by disk winds. Lee et al. ED, Fernan¬ 
dez & Metzger (SUES and Just et al EH have shown that 
5% — 25% of the matter in accretion disks formed in binary 
neutron star or black hole-neutron star mergers is eventually 
unbound, mostly due to viscous heating in the disk. In those 
2D simulations, and in the absence of a long-lived hypermas- 
sive neutron star, these outflows remained too neutron-rich 
to significantly affect the electromagnetic emission resulting 
from r-process nucleosynthesis in the ejecta (i.e. their nucle- 


FIG. 7. Temperature in the equatorial plane of the black hole at 
the end of the simulation, 20 ms after merger. Density contours are 
shown as thick black lines for po = 10 9 ’ 10,11 g/cm 3 . 



FIG. 8. Temperature in a plane orthogonal to the equatorial plane of 
the disk the at the end of the simulation, 20 ms after merger. Density 
contours are shown as thick black lines for po = 10 9,10,11 g/cm 3 . 


osynthesis output does not significantly differ from the out¬ 
put of the more massive dynamical ejecta) (52). However, 
the initial conditions for these simulations were taken to be 
Y e ~ 0.1 everywhere, at a time at which the disk is already 
wider and the composition evolution due to neutrino emission 
slower than in the simulations presented here. It is possible 
that a higher initial electron fraction in the outer regions of 
the disk can have a measurable impact on the properties of 
kilonovae. We discuss this in more detail in Sec. IIVDI Fi¬ 
nally, we consider the material in the disk winds. By the end 
of the simulation, these outflows have a high electron fraction, 
with Y e rsj 0.4. However, as discussed in Sec. |IVDl most of 
the material unbound in that region is ejected at earlier times, 
and is more neutron-rich. 
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FIG. 9. Electron fraction Y e in the equatorial plane of the black hole 
at the end of the simulation, 20 ms after merger. Density contours 
are shown as thick black lines for po — 10 9 ’ 10,11 g/cm 3 


in which the disk begins to deviate from hydrostatic equilib¬ 
rium, the Soldberg-Hoiland criterion is not strictly applicable. 
But the presence of convective regions is clearly observable at 
all times in the simulation. These instabilities cause the cre¬ 
ation of large scale eddies close to the outer edge of the disk, 
and are strongly correlated with outflows launched above the 
disk. Given how out of equilibrium the system is, picking a 
single instability responsible for the creation of the eddies or 
causally associating these eddies with the outflows is difficult. 
However, we note that regions in which there no longer are 
any interactions between the accretion disk and the tidal tail 
are devoid of both eddies and outflows. We can infer that the 
ejection of disk material is probably helped by the convec¬ 
tion of hotter fluid from the core of the disk to the top of the 
disk/tail interface. The high neutrino fluxes in that same re¬ 
gion probably play a role as well, especially in directing the 
outflows at late times and setting their composition. Compar¬ 
ing Fig. [6] and Fig. 12 it is clear that the disk outflows are, 
at the end of the simulation, aligned with the neutrino radia¬ 
tion. This is not the case at earlier times, when the outflows 
are stronger and closer to the spin axis of the black hole. 
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FIG. 10. Electron fraction Y e in a plane orthogonal to the equa¬ 
torial plane of the disk the at the end of the simulation, 20 ms af¬ 
ter merger. Density contours are shown as thick black lines for 
po = IQ 9,10,11 g/cm 3 . 


C. Stability and disk/tail interactions 

From the beginning of the simulation, strong asymmetries 
and shocks are present in the disk, and the system is well out 
of equilibrium. However, as the disk circularizes, it reaches 
a more stable configuration. By the end of the simulation, al¬ 
though asymmetries remain, fluid elements follow nearly cir¬ 
cular trajectories at the expected orbital frequency, up to radii 
r ^ 100 km. The core of the disk is convectively stable, and 
close to equilibrium. Around the disk/tail interface, however, 
this is not the case. First, the tail material is not circularized, 
which creates a shear layer at the outer edge of the disk. Ad¬ 
ditionally, according to the Soldberg-Hoiland criterion El, 
the disk is convectively unstable in both the vertical and ra¬ 
dial directions. We should note that as this is also the region 


D. Unbound material 

The ejection of unbound material by black hole-neutron star 
and neutron star-neutron star mergers is, from an astrophysical 
perspective, one of the most important consequences of these 
mergers. This is because these ejecta are one of the most likely 
locations for r-process nucleosynthesis to occur cuini Ei. 
Neutron star mergers were long thought to happen too late in 
the evolution of a galaxy to explain observations of r-process 
elements ED, but recent studies incorporating updated popu¬ 
lation synthesis predictions for neutron star binaries are more 
favorable to the merger process |58, t59|. Additionally, the 
radioactive decay of nuclei in the neutron-rich ejecta dur¬ 
ing r-process nucleosynthesis powers electromagnetic tran¬ 
sients potentially detectable in the optical or, more likely, in 
the infrared lfl4l [60l |6H . The luminosity, duration and peak 
frequency of these electromagnetic signals can significantly 
vary with the mass, composition, entropy and velocity of the 
ejecta El ED El. In particular, the composition will signif¬ 
icantly affect the results of the nucleosynthesis: low Y e ma¬ 
terial will produce mostly heavy elements (strong r-process), 
while high Y e material will produce more iron-peak elements 
(weak r-process). The transition occurs at Y e ~ 0.25 — 0.3 
for conditions typical of a viscously-driven wind in a post¬ 
merger accretion disk |[63l . However, nucleosynthesis results 
also depend on the entropy of the ejected material, and the 
exact dividing point for the lower entropy ejecta observed in 
our simulations is not, to our knowledge, known at this point. 
This difference in the products of r-process nucleosynthesis 
impacts the lightcurve of radioactively powered electromag¬ 
netic transients. Indeed, high-opacity lanthanides are pro¬ 
duces in the case of a strong r-process, causing the emission to 
be fainter, redder and longer lived than in the case of a weak 
r-process [14, 6JJ. 

In black hole-neutron star binaries, matter can be ejected 












12 


through different processes during and after the merger. First, 
if the neutron star is disrupted, a mass M e j ~ O.O1M 0 — 
0.1 M e (depending on the parameters of the binary) is typ¬ 
ically unbound during the tidal disruption of the neutron 
star MMM- The amount of ejected material can even 
be larger for rapidly spinning, low mass black holes [ 23, 46]. 
The tidally ejected material is very neutron rich (Y e < 0.1), 
cold (T < 1 MeV), confined close to the equatorial plane, and 
strongly asymmetric. 

After that, material can be ejected during the formation of 
the accretion disk and at the interface between the disk and 
the tidal tail. This is the main source of outflows observed in 
this simulation, and is discussed in more detail below. 

Third, disk winds may be triggered by magnetic effects 
and/or neutrino absorption. Newtonian simulations of bi¬ 
nary neutron star mergers using an energy-dependent leak¬ 
age scheme (with ad-hoc neutrino absorption) |[65l or energy- 
dependent flux-limited diffusion lf 66 l observed the formation 
of a neutrino-driven wind within ^100 ms of the merger, and 
neutrino driven winds were also observed in two-dimensional 
simulations of an accretion disk, with initial condition taken 
from a black hole-neutron star merger 1551 . In ESI, the ob¬ 
served outflows were more neutron-rich in the equatorial re¬ 
gion, with conditions favorable to a strong r-process, and less 
neutron-rich in the polar region, with the expectation of a 
weaker r-process (and bluer kilonovae). In ED, a wide range 
of electron fraction is observed in all directions, although the 
polar outflows remain less neutron rich. The disk generated 
in our black hole-neutron star mergers are likely, over longer 
timescales than those simulated here, to create winds similar 
to those observed in [[54 , 65 j. In those studies, about 1% of the 
disk mass was eventually ejected in the neutrino-driven wind. 
In black hole-neutron star mergers, this is generally much less 
material than in the dynamical ejecta. Accordingly, that ejecta 
is only interesting if its different composition has observable 
consequences. 

Finally, over longer timescales (seconds), and using an ide¬ 
alized initial disk profile, 2D Newtonian simulations have also 
shown that viscous heating can drive strong outflows in the 
disk H52j|67). The total ejected mass is about 5% — 25% 
of the mass of the disk, with some significant dependance 
on the spin of the black hole l53ll . Starting from idealized 
initial conditions (equilibrium torus with Y e ~ 0.1), these 
simulations find that most of the unbound material remains 
at electron fractions too low to avoid the production of lan¬ 
thanides - and they thus lead to an electromagnetic signal 
peaking in the infrared, and to strong r-process nucleosynthe¬ 
sis. For an initial black hole spin similar to our simulation 
(Xbh = 0.8), however, 0.1% — 1% of the mass of the disk is 
ejected in a less neutron-rich outflow, which does not produce 
any lanthanides and could lead to a ‘blue bump’ in the kilo- 
nova lightcurve ED. More recent results using more accurate 
neutrino methods also find large ejected masses, of the order 
of 20% — 25% of the mass of the disk [ 54]. 

In our simulation, we observe nearly-polar outflows during 
the circularization of the disk, mostly coming from the region 
in which the accretion disk and the tidal tail interact. The total 
mass ejected is only 3 x 1O _ 4 M 0 ~ 0.4%Mdi S k> a negligible 


amount compared to the mass ejected during the disruption 
of the neutron star (O.O 6 M 0 ). But because it is ejected early, 
and in a direction in which no material has been unbound so 
far, its effects cannot be neglected. First, material in the po¬ 
lar region can affect the formation and collimation of a jet, if 
the merger leads to a short gamma-ray burst. Additionally, be¬ 
cause this material is ejected before any disk wind, it could ob¬ 
scure the blue component of a kilonova if enough high-opacity 
lanthanides are formed during r-process nucleosynthesis and 
the ejected material obscures a significant fraction of the po¬ 
lar regions. The properties of the outflows observed in our 
simulation vary significantly in time. The front of the outflow 
has a relatively low electron fraction (Y e < 0 . 2 ) and specific 
entropy s ~ 30k b per baryon. At later times, the outflows 
become less neutron rich and colder. By the end of the simu¬ 
lation, we measure Y e > 0.3 and s ~ 15 ks- Overall, about 
15% of the ejected material has Y e < 0.2, and about 15% has 
Y e 0.3. Accordingly, we would expect most of the material 
unbound during disk formation to undergo strong r-process 
nucleosynthesis, and produce lanthanides. By the end of the 
simulation, the mass loss to these outflows has stabilized to 
about O.OO5M 0 /s. It is thus conceivable that after the end of 
our simulation, another ~ 1O - 4 M 0 of less neutron-rich ma¬ 
terial would be unbound through the same process. 

Putting all potential components of the ejecta together, we 
can get estimates of the various components of the ejecta for 
the specific binary parameters studied here. The quantitative 
results will of course vary with the initial parameters of the 
binary, but we expect the main qualitative features to be fairly 
robust within the most likely range of black hole and neutron 
star masses, at least as long as the black hole spin is large 
enough for the neutron star to be disrupted before merger. 
Here, we find: 

• About O.O 6 M 0 ejected asymmetrically in the equato¬ 
rial plane at the time of tidal disruption. This material 
is cold, neutron rich, and certainly leads to strong r- 
process nucleosynthesis and electromagnetic emission 
peaking in the infrared. 

• About 3 x 1O - 4 M 0 ejected at higher latitudes on a 
timescale of about 10 ms after merger, due to the circu¬ 
larization of the disk. This material is less neutron rich 
and hotter, with more uncertain nucleosynthesis output 
/ lightcurves (possibly depending on the initial param¬ 
eters of the binary). For the configuration considered 
here, however, we expect most of that material to un¬ 
dergo strong r-process nucleosynthesis. 

• Probably about 1O - 3 M 0 ejected in neutrino-driven 
winds on a timescale of about 100 ms after the merger, 
if our disk behaves as for the lower mass system studied 
in (65|. About half of that material should be ejected in 
the polar direction, with a higher Y e / entropy than for 
the material ejected in our simulation. The rest is closer 
to the equatorial plane, and more neutron rich. 

• Probably about 0.01 M 0 of material unbound due to 
viscous heating, on timescales of the order of 1 s [52j 
MWLl Most of that material is too neutron-rich to 
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avoid strong r-process nucleosynthesis, but a small frac¬ 
tion (rsj 10 _ 3 M 0 ) could be unbound in a hotter, less 
neutron rich front (63). 

We would thus end up with about 0.07M 0 of neutron-rich 
material ejected either in the equatorial regions or late in the 
polar regions. In the polar region, that material is preceded by 
ejecta with a much more complex composition. On the order 
of 10 _ 3 M q of less neutron-rich material is ejected, with the 
potential to produce different nucleosynthesis. If it is not ob¬ 
scured by neutron-rich material, it can also produce different 
kilonovae lightcurves (i.e. a blue bump). However, electro¬ 
magnetic radiation from that material could be significantly 
obscured by the lower Y e outflows observed in our simulation. 
Accordingly, the spatial distribution and exact mass of the 
outflows produced during disk formation probably should not 
be neglected when studying kilonova lightcurves from black 
hole-neutron star mergers. Considering that we only studied 
one configuration, and that the ejected mass observed here in 
the polar regions (3 x 10 _ 4 M o ) is barely large enough to mat¬ 
ter, it would be dangerous to draw overly generic conclusions 
from our results. The initial parameters of the binary, which 
significantly affect the distribution of matter between the disk 
and tidal tail and the conditions under which the disk form, 
are very likely to cause variations in the properties of these 
early disk outlows. But in the configuration studied here, the 
outflows cover a fairly wide part of the high-latitude regions, 
roughly 0 < 45° with 0 the angle with respect to the spin axis 
of the black hole. 


E. Neutrino Emission 



FIG. 11. Neutrino spheres r = 2/3 at the end of the simulation, 
20 ms after merger. We show the vertical plane y — ^bh. The three 
contours are, from the outside to the inside, the neutrinospheres for 
v e (outer solid line), v e (dashed line) and v x (inner solid line). The 
disk is colored according to its baryon density po- 
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With our Ml code, we can for the first time examine the 
spatial distribution of neutrinos in the first 20 ms following a 
black hole-neutron star merger. To understand the main prop¬ 
erties of the neutrino radiation, it is however useful to take 
a step back and look at a few of the quantities which could 
already be predicted using our leakage scheme. Indeed, the 
leakage scheme gives us a good estimate of the optical depth 
in the disk, a useful quantity to understand where the neutrinos 
effectively decouple from the fluid. For the electron neutrinos, 
although there are always a few hotter regions in which the op¬ 
tical depth is r Ve ^ 10 , the density averaged optical depth is 
only (r Ue ) ~ 2. It is even lower for the electron antineutrinos, 
with ( T ^e) ~ and heavy lepton neutrinos, with ( T v x ) ~ °- 5 - 
The neutrinospheres, defined as r = 2/3, are shown for a ver¬ 
tical cut at the end of the simulation in Fig. m We see that 
neutrinos are only trapped in the core of the disk, and mostly 
free streaming everywhere else. 

To illustrate the main properties of the neutrino radiation, 
we plot in Figs. [T2|[T3| the energy density and “normalized 
flux” F* orm = olF 1 /E — in vertical and horizontal slices 
of the disk, for the electron neutrinos. The normalized flux 
is chosen so that it represents an effective transport velocity 
for the neutrino energy. In the core of the disk, Fig. [13] shows 
that the neutrino energy is mostly transported with the fluid, 
as befits an optically thick region. Outside of the expected 
neutrinosphere, i.e. for r > 90 km, the neutrinos transition to 


FIG. 12. Energy density and normalized flux (aF l /E — f3 l ) of the 
electron neutrinos at the end of the simulation, 20 ms after merger. 
We show the vertical plane y — ?/bh . 


free streaming away from the disk. The energy density is max¬ 
imal close to the inner edge of the disk, in part due to higher 
temperatures and in part due to gravitational redshifting. 

The vertical slice on Fig. [12] shows a few additional fea¬ 
tures of interest. First, we note the random orientation of the 
fluxes along the speed axis of the black hole, up to a height 
\z\ < 8OM 0 . This is a known problematic feature of the Ml 
approximation, due to the convergence of neutrinos from all 
around the disk. Beyond this issue, we can also see that the 
emission at large radii is clearly asymmetric. We find signifi¬ 
cantly lower energy densities within about 20 ° of the equato¬ 
rial plane as well as in the (less reliable) polar regions. This is 
a purely geometric effect due to the projected shadow of the 
disk and the relativistic beaming of the neutrinos. 

To understand this beaming, we first have to remember that 
the velocities in the disk are mildly relativistic, with v ~ 0.5c. 
This means that the emission of neutrinos will be focused 
within a relatively large beam centered on the direction of mo¬ 
tion of the fluid, which nearly follows a circular orbit around 
the black hole. Only a small fraction of the neutrinos are emit¬ 
ted in the vertical direction. Additionally, most of the neutri- 
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FIG. 13. Energy density and normalized flux (aF l /E — /3 Z ) of the 
electron neutrinos at the end of the simulation (20 ms after merger), 
in the equatorial plane of the disk. 



The radiation fields of the other neutrino species are qual¬ 
itatively similar, except for the fact that the neutrino sphere is 
located deeper inside the disk. We find that most of the energy 
emitted in neutrinos goes into electron antineutrinos, at least 
at early times. The luminosity in v e is fairly constant during 
the evolution, at L v& ^ 5 —8 x 10 52 ergs/s. But the luminosity 
in v e peaks at L„ e ~ 5 x 10 53 ergs/s during the rapid evolu¬ 
tion of the electron fraction from its very neutron-rich initial 
value to the equilibrium value in the disk. Within 5 ms, it de¬ 
creases to Ly e ~ 10 53 ergs/s, and by the end of the simulation 
it becomes lower than the electron neutrino luminosity, with 
L[ e = 7 x 10 52 ergs/s and L^ e = 5 x 10 52 ergs/s. Finally, 
there is an equally brief burst of heavy-lepton neutrinos, with 
Ljy x ~ 3 x 10 52 ergs/s for each species, due to the existence of 
hot spots early in the disk formation. That luminosity rapidly 
decreases as the temperature becomes more homogeneous. 
Within 3 ms, we measure L Vx ~ 5 x 10 51 ergs/s for each 
species, and by the end of the simulation, L? = 10 51 ergs/s 
per species. 

The properties of the neutrino spectrum are not directly 
measurable within our gray formalism. Nonetheless, it is pos¬ 
sible to get reasonable estimates from either the predictions of 
the leakage scheme or the properties of the fluid on the neu- 
trinosphere. The leakage scheme predicts average energies of 
(e Ue ) = 11 MeV — 13 MeV, (e„ e ) = 14MeV - 15 MeV, and 
\e Vx ) = 16 MeV —18 MeV during the last 10 ms of evolution. 
All of the energies peak ^ 4 MeV higher at the beginning of 
the simulation, when hot spots are present. Neglecting cor¬ 
rections due to the finite chemical potential of neutrinos in 
the emitting regions and assuming a redshift factor of 2 be¬ 
tween the emitting region and the observer, this would cor¬ 
respond to fluid temperatures of roughly 6 MeV, 7 MeV, and 
8 MeV. Considering the temperature distribution observed in 
the disk, the increasing temperature as the neutrino sphere re¬ 
cedes deeper into the disk, and the higher emissivity of high 
temperature points, this appears fairly reasonable. 


FIG. 14. Energy density and normalized fluxes of the electron neu¬ 
trinos on the r Ue =0.1 surface. F. Neutrino viscosity and the growth of magnetic instabilities 


nos decouple from the matter around the location of the neu- 
trinosphere, and then free stream away from the disk. Accord¬ 
ingly, if the inner regions of the disk are brighter than its outer 
regions, as is the case here, the disk casts a shadow along the 
equatorial plane. This is more clearly illustrated by Fig. |T4| 
which shows the energy density and fluxes for the electron 
neutrinos on the surface r Ue =0.1. Most of the energy comes 
from the inner disk, and neutrinos cannot easily escape along 
the equatorial plane. The shadow along the equatorial plane 
observed in Fig. [12] exactly matches the thickness of the disk. 
Relativistic beaming, on the other hand, causes neutrinos to 
be preferentially emitted nearly tangent to the disk. Hence the 
neutrino flux is larger just outside of the shadow than at high 
latitudes. The equatorial shadow is of course not perfect, as 
some neutrinos are emitted from the outer edge of the neutri- 
nosphere. But this is a relatively small fraction of the total 
neutrino luminosity, as can be seen in Fig. [14] 


The magnetorotational instability (MRI) is expected to play 
a crucial role in the growth of magnetic fields in post-merger 
accretion disk, and may also be important for the genera¬ 
tion of jets after a compact binary merger. However, recent 
work has shown that the transport of momentum by neutri¬ 
nos can significantly affect the growth timescale and wave¬ 
length of the MRI ( 68 ] [69). In the context of protoneutron 
stars, Guilet et al. [69) showed that if the wavelength of the 
fastest growing mode of the MRI in the absence of neutrinos 
Amri is larger than the neutrino mean free path A^, an effec¬ 
tive viscosity from the transport of neutrinos causes the MRI 
to grow slower than expected for small magnetic fields, but 
at a fastest-growing wavelength independent of the magnetic 
field strength (instead of the wavelength decreasing with the 
magnetic field strength). If instead Amri < A^, the neutri¬ 
nos can act as a drag force on the magnetized fluid. For large 
neutrino energy densities, this slows the growth timescale of 
the fastest growing mode of the MRI, and slightly increases 
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its wavelength. In accretion disks, and when neutrinos can 
be modeled through an effective viscosity, a similar increase 
of the growth timescale of the MRI has been measured [ 68]. 
By applying the model of Guilet et al. l69l to our accretion 
disk, we can obtain a first estimate of the expected effect of 
neutrinos on the growth of the MRI. 

First, we need to determine the critical magnetic field B c at 
which A M ri = K, using A M ri ~ 27 rb/>Jb 2 + p 0 h and A^ ~ 
1 /+ k s ) (where b is the strength of the magnetic field ob¬ 
served by an observer comoving with the fluid). At the high¬ 
est density points in the disk, we get B c ~ 10 13 G, while B c 
smoothly increases with decreasing density, to B c ~ 10 15 G 
for po ~ 10 10 g/cm 3 . These values are larger than the initial 
magnetic field in most merging neutron stars, but smaller than 
the expected saturation amplitude of the MRI, thus indicat¬ 
ing that both the viscous and neutrino drag regimes might be 
relevant during the evolution of a post-merger accretion disk. 

As far as the neutrino drag is concerned (e.g. for B < B c ), 
its effects on the growth of the MRI should however be negli¬ 
gible. Guilet et al. l69l find that the importance of neutrinos 
on the growth of the MRI in this regime is determined by the 
value of the dimensionless parameter 


r _ 2(ft a -f- n s )E v 

Q lbpoQ 

For T/fi > 1, neutrinos affect the growth of the MRI. How¬ 
ever, in our disk, we find T/O < 0.01 everywhere. Thus the 
growth of the MRI should remain unaffected at low magnetic 
field strengths. 

In the viscous regime (B > B c ), the importance of neutrino 
effects is determined by the value of the Elasser number e = 
v\/(afl) where va is the Alfven speed and 


2 E v 

15po(^a + Ks) 


(35) 


is the effective viscosity due to neutrinos. Viscosity affects the 
growth of the MRI for E u < 1. For B ~ B c , and the condi¬ 
tions observed in our simulation, we would get e ~ 1 (and the 
Elasser number then grows as B 2 ). Accordingly, within the 
simple model used here, neutrinos could plausibly affect the 
growth of the MRI for B ~ B c . But for e ^ 1, these effects 
are mild, and neutrinos are not expected to have any effect for 
B < B c or B > B c . 


G. Impact of the neutrino treatment 

Having performed simulations with both a leakage scheme 
and the Ml formalism, we can now obtain better estimates of 
the impact that an approximate treatment of the neutrinos has 
on the evolution of a post-merger accretion disk. 


properties of the high density regions of the disk: the sim¬ 
ulations with neutrino leakage capture the formation of the 
disk, its early protonization and later re-neutronization, and 
the global temperature evolution due to the initial expansion 
of the disk, neutrino cooling and shock heating. Lacking 
neutrino absorption, the leakage simulation however tends to 
maintain larger temperature differences between neighboring 
regions of the disk. It underestimates the timescale necessary 
for the disk to cool down, and the magnitude of its protoniza¬ 
tion. More quantitatively, the average electron fraction in the 
leakage simulation is lower by A Y e ~ 0.05, and the average 
temperature is lower by AT = (0.5—0.8) MeV. Additionally, 
without neutrino absorption the evolution of the fluid compo¬ 
sition in low-density regions is entirely unreliable. The com¬ 
position of the tidal tail material falling back onto the black 
hole and of the outflows are thus radically different. 

There are also significant differences in the neutrino lumi¬ 
nosities, for the electron neutrinos by a factor of two about 
10 ms after merger (5 ms after the beginning of the simula¬ 
tion), and by about 30% by the end of the simulation. For 
heavy-type neutrinos, the difference is typically a factor of 
2 — 4, presumably because of the steeper dependence of their 
emissivities on the temperature of the disk. We show the lu¬ 
minosities for the leakage scheme, the Ml scheme, and the 
predictions of the leakage scheme for the value of the hydro- 
dynamic variables obtained when evolving the system with 
the Ml scheme on Fig. 15 The first thing to note is that, 
since the leakage scheme measures the instantaneous energy 
loss at each point of the domain while the luminosities in the 
Ml scheme are measured through the neutrino fluxes at the 
outer boundary of the computational domain, there is natu¬ 
rally a time shift between the two predictions. Even taking 
that time shift into account, however, it is quite clear that the 
electron neutrino luminosity is larger in the leakage scheme, 
while the heavy-lepton neutrino luminosity is larger in the 
Ml scheme. The electron antineutrino luminosities cannot be 
distinguished within the uncertainties due to the time shift in 
the measurements. Looking at the predictions of the leakage 
scheme within the Ml simulation allows us to differentiate 
between the error due to the instantaneous estimate of the lu¬ 
minosity, and the error due to the diverging evolution of the 
hydrodynamics variables. We see that by the end of the simu¬ 
lation, both sources of error are important. 

Overall, these results appear consistent with the expected 
limitations of a leakage scheme. We also performed leakage 
simulations with two different methods to compute the opac¬ 
ities, to check whether the differences observed between the 
leakage and Ml schemes could be due to differences in the es¬ 
timated neutrino opacities. But the differences between these 
two simulations were well below the estimated errors in the 
leakage scheme. 


1. Comparison with a leakage scheme 2 . Impact of the gray approximation 


We can first look at differences between the leakage simula¬ 
tions and the Ml simulations. Not surprisingly, the two meth¬ 
ods provide reasonable qualitative agreement for the global 


By looking at Ml simulations using different energy¬ 
averaging approximations, we can also attempt to estimate er¬ 
rors related to the use of a gray scheme. We use two different 
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FIG. 15. Luminosities for the different neutrino species as a function 
of time. Solid curves show the Ml results. Dashed curves show re¬ 
sults using a leakage scheme. Thin dot-dashed curves show the pre¬ 
diction of the leakage scheme, but for the value of the hydrodynamic 
variables obtained by evolving the disk with the Ml scheme. Heavy- 
lepton neutrino luminosities v x are for all 4 species combined. 


gray approximations which, in our test of neutrino-matter in¬ 
teractions in a post-bounce supernova profile, bracketed the 
solution obtained with an energy-dependent code (see Ap¬ 
pendix [E]). Clearly, as an error estimate, this is a poor sub¬ 
stitute for a comparison with an energy-dependent radiation 
transport code. Unfortunately, such a simulation remains too 
costly to perform with the SpEC code at this point. We find 
that differences in the average properties of the disk, although 
measurable, are much smaller than the differences between 
the Ml simulations and the leakage simulations. The average 
electron fraction varies by A Y e < 0.01, and the average tem¬ 
perature by AT <0.1 MeV. The total neutrino luminosities 
agree, for all species, within ~ 20%. 

We also performed a simulation in which the gray opaci¬ 
ties were computed assuming that the neutrinos obey a Fermi- 
Dirac distribution with temperature T equal to the fluid tem¬ 
perature and equilibrium chemical potential a method 
which significantly underestimates neutrino opacities: low 
density regions can have T ~ 1 — 2 MeV while neutrinos 
are mostly emitted in regions with T > 5 MeV. The result¬ 
ing differences are significantly smaller than in the test of the 
neutrino-matter coupling presented at the end of the appen¬ 
dices. In fact, in the core of the disk, the errors are similar to 
the differences between the two Ml simulations using differ¬ 
ent prescription for the emissivities. Where the two methods 
differ, unsurprisingly, is in the composition of the outflows: 
when assuming an equilibrium distribution at the fluid tem¬ 
perature for the neutrinos, the electron fraction of the outflows 
is consistently AV e = 0.05 lower than when using the more 
realistic neutrino energies computed from the temperature in 
the emitting regions. This confirms our assumption that the 
electron fraction in the outflows is significantly affected by 
neutrino absorption. 


3. Limitations of the Ml closure 

Finally, when analyzing our results, it is worth noting once 
more the limitations of the Ml closure. Crossing beams, caus¬ 
tics, and strongly focused beams are known to be problematic 
when using the Ml closure. Some examples of these issues 
can be found in the tests presented in Appendix [E] Accord¬ 
ingly, the results of our simulations cannot be trusted close to 
the spin axis of the black hole, where neutrinos from all re¬ 
gions of the disk cross. Because of the relativistic beaming of 
the neutrinos, the energy density in that region is already low, 
but unphysical radiation shocks cause an even larger decrease. 
Unfortunately, the neutrino radiation in the polar region can 
have important consequences on the post-merger evolution of 
the system. Indeed, neutrino-antineutrino annihilations into 
electron-positron pairs deposit energy in this low-density re¬ 
gion, and affect the matter density there. The formation of 
a relativistic jet, desirable if black hole-neutron star mergers 
are to produce short gamma-ray bursts, is quite sensitive to 
baryon-loading in the polar regions, but on the other hand 
could be helped by the energy deposition. The exact impact of 
the neutrinos in that context remains an open question, which 
can only be answered with an improved treatment of the neu¬ 
trinos, and a better understanding of the jet forming mecha¬ 
nism. 


V. CONCLUSIONS 

We present a new module of the SpEC code ED, allow¬ 
ing us to study the effects of neutrinos within a fully general 
relativistic - hydrodynamics code. The neutrinos are modeled 
using the Ml formalism, in which the first two moments of the 
neutrino distribution function are evolved (energy and fluxes). 
Although the formalism can in theory be energy-dependent, 
we limit ourselves to an energy-integrated version of the code, 
due to the high cost of energy-dependent simulations. We of¬ 
fer here a detailed description of the implementation of the Ml 
algorithm in SpEC, as well as a series of tests assessing our 
ability to study the evolution of neutrinos in flat and curved 
spacetimes, and their interaction with matter. 

We also discuss the first simulation with a general relativis¬ 
tic code of the evolution of an accretion disk produced by a 
typical black hole-neutron star merger. We use as initial con¬ 
ditions a snapshot of an existing SpEC simulation using a sim¬ 
pler treatment of the neutrinos (leakage scheme) t24ll . right as 
the accretion disk begins to form and neutrino effects become 
important. This provides us with realistic initial conditions for 
the accretion disk, at least compared to the more commonly 
used equilibrium tori. As we neglect magnetic fields and all 
material falling back on the disk on timescales longer than 
20 ms, we limit ourselves to a relatively short evolution, up to 
20 ms after merger. We find that the evolution of the forming 
accretion disk is initially dominated by the circularization of 
the disk material. The disk expands and contracts in a cycle of 
about 6 ms. These oscillations are the main driver of the evo¬ 
lution of the average disk properties at early times. As the disk 
circularizes, however, these oscillations are rapidly damped. 
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Neutrinos cause a global cooling of the disk, but on a rela¬ 
tively long timescale of order 50 ms. On the timescale of this 
simulation, neutrinos have however other important effects. 
First, they drive the composition evolution of the core of the 
disk. The disk rapidly protonizes, reaching an average elec¬ 
tron fraction Y e ~ 0.2 (a larger value than predicted by the 
leakage scheme [ 241). The electron fraction then more slowly 
decreases as the disk starts to cool down. Second, neutrinos 
cause changes in the electron fraction of the outer parts of the 
disk, and of the disk outflows. Although initially neutron rich, 
by the end of the simulations these regions have Y e > 0.3 
(Y e rsj 0.4 in the outflows). These low density regions above 
the disk would be the first to be unbound by disk winds and 

viscously-driven outflows E2H33II32, and their higher initial 

electron fraction could thus affect the composition of late time 
outflows, their nucleosynthesis output, and the light curve of 
associated electromagnetic transients. Third, energy transport 
by the neutrinos homogenize the temperature distribution in 
the disk, with cold regions absorbing higher energy neutrinos 
and hot regions radiating much more strongly. Finally, neu¬ 
trinos will deposit energy and create electron-positron pairs 
in the low-density polar regions. This could plausibly affect, 
positively or negatively, the formation of relativistic jets in 
that region. Unfortunately, neutrino-antineutrino annihilation 
in low density regions is not modeled within our formalism. 
As opposed to what was observed in simulations of protoneu¬ 
tron stars l69l . we also estimate that the impact of neutrino 
transport on the growth of the magnetorotational instability is 
likely to be small to inexistant in post-merger accretion disks. 

We also find that, because the disk is very compact, and 
forms around a rapidly rotating black hole (xbh = 0.87), rel¬ 
ativistic effects are significant. Most of the radiation comes 
from a region with mildly relativistic velocities (v ~ 0.5c). 
Accordingly, the neutrinos are beamed tangentially to the 
disk, causing strong anisotropies in the neutrino luminosity. 
Light bending and gravitational redshift also naturally affect 
the luminosity and neutrino spectrum. The luminosity is low 
in the equatorial plane, due to the disk shadow, and in the po¬ 
lar regions, due to beaming (although the evolution of polar 
regions with the Ml closure is unreliable). We note that even 
though a black hole spin %bh = 0.87 might seem large, this 
should be fairly typical of a black hole-neutron star merger 
in which an accretion disk forms. Indeed, for the most likely 
black hole masses, slowly spinning black holes cannot disrupt 
the neutron star before it plunges into the black hole sa. 

A number of simulations using approximate treatments of 
gravity have also considered the impact of neutrinos on disk 
evolutions, with methods generally more advanced than con¬ 
temporary general relativistic simulations (e.g. El EH m 
l70l I7T]), and/or capable to evolve the system over longer 
timescales (e.g. [51 4531). However, direct comparisons with 
our results are difficult. This is in part because of the impor¬ 
tance of relativistic effects for accretion disks around rapidly 
spinning black holes, and also because simulations which do 
not include general relativistic effects typically start from ei¬ 
ther an idealized equilibrium torus, or from the result of a 
merger in which the neutron star was disrupted by tidal ef¬ 
fects modeled by a pseudo-Newtonian potential, whose real¬ 


ism close to a rapidly spinning black hole is difficult to assess. 
Our simulation shows that general relativistic effects signifi¬ 
cantly impact the neutrino radiation and the disk formation, 
and the forming accretion disk is more compact and hotter 
than in non-relativistic studies. An important application of 
our results would in fact be to provide better initial conditions 
for long-term disk evolutions, or for nucleosynthesis studies 
requiring a detailed knowledge of the disk structure and com¬ 
position (e.g. l72l ). 

Finally, we note that the joint effects of shocks during the 
disk circularization, instabilities at the disk/tail interface, and 
neutrino absorption unbinds a small amount of material in the 
polar regions (~ 3x 10 _4 M o ). This might seem negligible 
compared to the material ejected dynamically in the equatorial 
plane during the disruption of the neutron star (~ 0.06M o ). 
However, this ejecta could be important because it is unbound 
in a direction in which no material has been ejected so far, 
and could thus impact the formation of a relativistic jet. Addi¬ 
tionally, over longer timescales, we expect neutrino-powered 
winds to become active and eject material in the polar regions, 
maybe of the order of 1% of the mass of the disk. Later 
on, viscously-driven outflows could eject a more significant 
amount of material (probably 5% — 25% of the disk). But the 
outflows observed here will be the outermost layer of ejected 
material in the polar region. Their opacity might affect the 
properties of electromagnetic transients for observers in direc¬ 
tions close to the spin axis of the black hole. From our simu¬ 
lations, it appears that most of the matter in these outflows is 
too neutron rich and cold to avoid strong r-process nucleosyn¬ 
thesis and the formation of high-opacity lanthanides, and thus 
that these early outflows could obscure later disk winds. How¬ 
ever, we caution that the ejected mass and geometry of the 
outflows are likely to depend on the parameter of the binary, 
and thus that it would be dangerous to draw overly generic 
conclusions from a single initial configuration. 

With these results demonstrating the ability of the SpEC 
code to evolve neutrinos within the moment formalism, we 
are now in a position to improve on the short simulations pre¬ 
sented here. The gray scheme used in this work, which would 
be inadequate for the study of core collapse supernovae, ap¬ 
pears at this point sufficient for the study of neutron star merg¬ 
ers. 
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Appendix A: Closure relation 


The choice of a closure P^ ( E , F^) is the main approxima¬ 
tion used in the moment formalism. In this work, we use the 
Ml closure, which relies on an interpolation between the ex¬ 
pected closure relation in the optically thin and optically thick 
regimes. Although we evolve E and Fi, the closure relation 
is generally easier to express as a function of the fluid frame 
quantities J and H Following (34), we define the parameter 
c as 


H«H a 

J 2 5 


(Al) 


such that ( ~ 0 in the optically thick limit and ( ~ 1 in the 
optically thin limit. The closure is then 


pij = 


3x(C) - 1 pij 

o ^thin 


3(1 -X(0) 

2 


pij 

r thick ' 


By default, we use the Minerbo closure l36l 


(A2) 


x(C) 



6 — 2C + 6 C 2 


(A3) 
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In the optically thin limit, we then use 


du 

^thin 


F3 


E 


(A4) 


Although this expression is exact in the limit E 2 = F l Fi , 
Shibata et al. El have shown that it does not respect causality 
when E 2 > F l Fi. However, proposed alternatives have more 
serious limitations El, and P t ^ in only matters when E 2 ^ 
P'Pi. 

For the optically thick limit, we choose 


where e v is the neutrino energy, fi v the neutrino chemical po¬ 
tential, ks the Boltzmann constant, and T the temperature of 
the fluid. 

We now have to choose the value of the chemical potentials 
/^(po, T,Y e ). The value of the chemical potentials for neu¬ 
trinos in equilibrium with the fluid, p® q , is taken directly from 
the equation of state table. We consider two different choices 
for ji u : either we set ji u = everywhere, or we make the 
same choice as in many leakage codes, 

= l-e~ r ") , (B2) 


q^ v — 

D thick ~ n 


^thick 


(A5) 


which is equivalent to 


Jthick = 2W2 + 1 [(2^ 2 - l ) E - 2W 2 F i Vi] , (A 6 ) 

7 ^< i Ck = P + [( 4 ^ 2 * + ')*** ~ * w2p ] • 


P t hick can ^ en computed from Eq. ( 14}. 


In practice, because P l ° is a function of £, which itself de¬ 
pends on J, H a , which are themselves functions of P 2 - 7 , the 
equations that we just described can only be solved through 
the use of a root-finding algorithm. We thus define 


P(C;P,Pi) = 


CJ 2 -H«H a 
E 2 


(A7) 


and solve for P(£; E, Fi) =0 using a Newton-Raphson algo¬ 
rithm (with ( initialized at its last computed value at the given 
point). In Eq. m J and H fJ are computed from P, F, and 
Pij , where P^ is now itself a known function of E , Fi and 


Appendix B: Energy integrated source terms 

To evolve the energy-integrated moments E and Fi of the 
neutrino distribution function , we need to define the emis- 
sivity rj and opacities (/s a ,ft s ) of the fluid. In theory, these 
source terms are functions of f^y However, as we only 
evolve energy-integrated moments, we do not know the neu¬ 
trino spectrum, and have instead to rely on estimates of . 
There are multiple ways to do so, and we detail the various 
choices that we have implemented below. We note that a num¬ 
ber of those choices rely on information obtained from a sim¬ 
pler leakage scheme for the treatment of neutrino cooling |[23l 
(e.g. average energy of neutrinos, optical depth), so that even 
when we evolve the neutrino radiation using the moment for¬ 
malism, we still leave the leakage scheme turned on — but 
without coupling it to the fluid equations. 

We can first compute the source terms assuming that the 
neutrinos obey a Fermi-Dirac distribution with temperature T 
and chemical potential fi v . In that case, we have 

^ 1 • exp|(<=„ - //„)/(A-, ; '/')| 


which is chosen so that /i v 0 in the optically thin region. 
We will refer to these two choices as ‘equilibrium’ and ‘leak¬ 
age’ chemical potentials. The choice of fi v in the low optical 
depth regions does not matter much in practice when using 
our default method for the computation of the source terms 
(described below), but can have an effect in some of the alter¬ 
native schemes that we tested. 

The energy integrated source terms are then computed as 
in our leakage code l23l (which is itself based on the GRID 
code 139], and previous work by Ruffert et al. E3, Rosswog 
and Liebendorfer ESI and Burrows et al. El): we compute 
the absorption opacity due to electron and positron capture, 
the free streaming emission due to e + e - pair annihilation, 
plasmon decay and nucleon-nucleon Bremsstrahlung, and the 
scattering opacities due to elastic scattering on nucleons and 
heavy nuclei. We note however that, for the emissivity in opti¬ 
cally thick regions, we use the free streaming emissivity, while 
the leakage scheme replaces the emissivity by the diffusion 
rate in those regions. Indeed, diffusion through the optically 
thick regions is handled naturally when using the moment for¬ 
malism. 

In order to guarantee that the neutrinos are in equilibrium 
with the fluid in the optically thick regions, we then com¬ 
pute the free streaming emission due to charge current re¬ 
actions and the absorption due to pair processes through an 
energy-integrated version of Kirchhoff’s law. At a given en¬ 
ergy Kirchhoff’s law gives us a relation between the emis¬ 
sivity rj (e^), the absorption opacity n a {e y ), and the equilib¬ 
rium spectrum of the neutrino radiation B v {e v >0 

V(ev) = K a (e v )B u (e v ) . (B3) 


The absorption opacity entering the evolution equations for E 
and Fi, on the other hand, is the energy averaged 


Jq (e^)/(e^)de^ 

r 


(B4) 


where I(e u ) is the specific intensity of the neutrino radiation. 
In our code, when computing K a from r] (for pair processes), 
or r] from K a (for charged-current reactions), we assume that 


5 Note that Eq. |B3| is always true for charged-current reactions, but only in 

the equilibrium limit for pair processes. 
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I u = B v . Then, we have 

POO POO 

V= v(tv)de v = / K,a(e„)B u (e u ) 

Jo Jo 

oo 

B v {e v )de v . (B5) 

As B v is a known function of the fluid properties, this equa¬ 
tion can easily be enforced. In optically thick regions, this 
prescription will be accurate, as it maintains the equilibrium 
between the fluid and the neutrino radiation. In optically thin 
regions, on the other hand, the neutrino radiation can be far 
from equilibrium. Computing charged-current interactions 
from the energy-integrated Kirchhoff’s law assuming an equi¬ 
librium spectrum can affect the relative emission of electron 
neutrinos and antineutrinos, and thus the evolution of the elec¬ 
tron fraction in low-density regions. We thus consider an alter¬ 
native to the application of Kirchhoff’s law when computing 
the emissivity of charged-current reactions: we can smoothly 
interpolate between the value of the emissivity t\k predicted 
by the application of Kirchhoff’s law and the emissivity 77f ree 
predicted by Ruffert et al. ED for free emission in an optically 
thin region, using 

V = VKf{T v ) + 7?free[l ~ /(r„)] , (B6) 



with 


f( T v) = 1 if r v > 2 

= 0 if tv < 2/3 

= Ti/ otherwise . (B7) 

Note that even when using this corrected emissivity, the opac¬ 
ity due to charged-current reactions is not modified. Accord¬ 
ingly, this corrected emissivity no longer satisfies Eq. ( |B5| ). 
Although apparently more ad-hoc than the previous method, 
this prescription gives the best results in our test of neutrino- 
matter interactions (see Appendix 0. This is probably be¬ 
cause the energy integrated emissivity ? 7 f ree was specifically 
meant to approximate the emissivity in optically thin regions. 
We note that the deviations from Kirchhoff’s law for the en¬ 
ergy integrated K a and 77f ree come from the approximations 
made in computing energy-averaged Pauli blocking factors, as 
well as from neglecting the electron rest mass and the differ¬ 
ence between the neutron and proton masses when integrating 
over neutrino energies. 

Finally, we can go one step further in improving our evo¬ 
lution scheme by noting that for most reactions considered 
here, including the dominant charged-current reactions, the 
absorption and scattering cross-sections scale like e^. This 
means that in low-temperature regions, where the fluid is 
much colder than the neutrinos (which, presumably, are emit¬ 
ted in higher temperature regions), we largely underestimate 
K a , s . We can thus apply the correction 


max • 


r'v,leak/ 

fluid/ 


(B8) 


where (e^ leak ) is the average of predicted by the leakage 
scheme (taken over the entire grid) while 


fluid 


-^5 ijlv) rp 2 


(B9) 


is the average of for neutrinos in equilibrium with the 
fluid. The average is weighted by the energy density of neu¬ 
trinos and not the number density of neutrinos, as we are in¬ 


terested in the energy-averaged K a ,s given by Eq. (B4). In 
our tests (see Appendix [EJ, we find that using this corrected 
value of K, a ,s significantly improves the agreement between 
the gray scheme used in this work and energy-dependent radi¬ 
ation transport. 

We thus leave open three choices when computing the 
source terms. The chemical potential can be obtained from 
its equilibrium value or using the ‘leakage’ prescription. The 
emissivity of charged-current reactions in the optically thin 
limit can be obtained from the opacities of the fluid using the 
energy-integrated version of Kirchhoff’s law or from a direct 
estimate of the emissivity. And finally, the opacities in low 
temperature regions can be obtained by assuming equilibrium 
between the neutrinos and the fluid, or by applying a correc¬ 
tion for the higher energy of the neutrinos. Our default algo¬ 
rithm is to apply the energy correction to to compute rj 
in optically thin regions from a direct estimate of the emissiv¬ 
ity, and to use the equilibrium chemical potential everywhere 
(although that latest choice has nearly no impact when com¬ 
bined with our prescription for rf). We should however note 
that, even though these choices lead to significant differences 
in our test based on the evolution of the post-bounce config¬ 
uration of a core collapse simulation, they are in much better 
agreement in the case of binary mergers. 


Appendix C: Computation of the fluxes at cell faces 

When computing the fluxes at cell faces for the evolution 
of E and F z , we first use shock-capturing methods to estimate 
the value of ( E,Fi/E) at cell faces from their value at cell 
centers. For the simulations presented in this paper, we use 
the 2nd order minmod reconstruction (see e.g. Elim al¬ 
though higher order reconstruction methods are available in 
SpEC (the fluid variables at cell faces, for example, are re¬ 
constructed using the WEN05 algorithm [77, 78]). These re¬ 
construction methods are used with both left- and right-biased 
stencils, leaving us with two reconstructed values Ul and Ur 
for the evolved variable U = (E,Fi) and FI and Fr for the 
fluxes F l = (aF l — /3 l E,aPj — /3 l Fj). The fluxes at cell 
interfaces are then approximated by the HLL formula 


c+Fl + c-Fr — c+C-(U R — Ul) 
c + + c_ 


(Cl) 


where c + and c_ are the absolute values of the largest right- 
and left-going characteristic speeds of the evolution system 
(or zero if there is no left/right going characteristic speeds). 
For the characteristic speeds Ai 5 2 at a cell interface along di- 
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rection d, we use 


v d 

p = a w 


r = sJFFUw 2 + 1) - 2W 2 p 2 , 

A - & n \ Fd \ 

1>thin “ ^ y/FWi ’ 

Ai,thick = min (-/3 d + -p d + p ) 


A 


■2, thin 


2W 2 + 1 
-/3 d + a F\ 


VFWi 

7 2 
UF 

, 3(1 -x(0) 

i z,thin \ 0 ' 


A 2 ,thick = min (-/3 d + + P) 

Aj = ~ 1 A 


(C2) 

(C3) 

(C4) 

(C5) 

(C6) 

(C7) 




where the characteristic speeds in the optically thick and thin 
limits are taken from l34l and the interpolation between the 
two regimes uses the same formula as for the closure. 

This prescription works well in optically thin regions. In 
the optically thick limit, however, these fluxes do not properly 
reproduce the diffusion rate of the neutrinos through the fluid. 
We thus apply corrections to the flux F f of E 179 1 


Fe, corr — &Fe + (1 — cl)Fe ,asym i (C9) 


with 


a = tanh - 


K,Ax d 

i+l/2 = \!(^a F ^s)i{^a 4" ^s)z+l 5 


(CIO) 

ceil) 


and where half-integer indices refer to values of the opacities 
at cell faces while integer indices refer to value of the opac¬ 
ities at cell centers. Here, d is the direction in which we are 

reconstructing, Ax d = yjgdd(Ax d rid ) 2 is the proper distance 

between two grid points along that direction, and Ax d rid the 
coordinate grid spacing along that direction. The asymptotic 
flux in the fluid rest frame, which corresponds to the flux in 
the diffusion limit, is 03 


— 0 ^a^thick 


(Cl 2) 


Using this and equation |A5| in equation [13| gives the asymp¬ 
totic observer frame flux 


W 2 av d Jt h i ck - /3 d 7) 

/ di 1 d i\ /—^^thick o\ 

-a —(7 + vv)y/i ■ (C13) 

Numerically, this term can be fairly complex to evaluate. It 
requires derivatives of the neutrino energy density in the fluid 
frame along all coordinate directions, estimated at cell faces. 
In practice, we compute the first and second term of Fe, asym 
separately. For the second term, which models the diffusion of 
neutrinos, all quantities are estimated by averaging the values 


Fe, asym — \fl 


at n eighb oring cell centers, except for ft, which is given by 
Eq. (Cll), and dJ/dx d , for which we use 


/ dJ \ _ JiE 1 Ji 

\dx d ) i+ 1/2 “ Ax d 


(C14) 


For the first term, which represents the advection of neutrinos 
with the fluid, we make separate estimates from the left and 
right states (Ul, Ur). For both states, we also compute the 
advection speed 


Cadv = -7 + 4a 2VP + U ' ( C15 ) 

If both speeds are positive, we use the value computed from 
Ur. If both are negative, we use the value computed from Ur. 
If their signs disagree, the advection term is set to zero. 

A simpler correction is also applied to the flux Fp of Fi, 
following |80l 


Ff,cott = A Fe + (1 - A ) 


2\Ff,R + Fp,L 


with 


A = 


1 

RAx d 


(Cl 6) 


(C17) 


Appendix D: Implicit time stepping 


The collisional source terms in equations ( I6p7 ) can be 
very stiff in optically thick regions. If we were to treat those 
source terms explicitly, we would need to use prohibitively 
small time steps. Accordingly, we will split the evolution 
equations into implicit and explicit terms. The variables 
(E n+1 , F iin+1 ) evolved from (E n , Fi^ n ) by taking a time step 
dt are given by 


+ dj{aFi-PE n ) (Dl) 

= a(Pj°K ij -F’d j ]na-SZ + 1 n a ) , 

+ dj(aPR - ?F itn ) - a~S« +llia (D2) 
= {-Endiot + F k , n di(3 k + ^P^dajk) . 


Solving this exactly would require the use of a four¬ 
dimensional non-linear root-finder algorithm to get 
(E n +i, F n +i). To limit the cost of each time step, we 
instead use a linear approximation to 5“ +1 , 


~Sn + i = + A a E n + 1 + B^F^ n+ 1 , (D3) 

where the coefficients A a and B are computed assuming 
that the closure factor £ and the angle between the neutrino 
flux Fi and the fluid velocity v l are constant. The system then 
reduces to four linear equations for (E n +i, i^ jn +i), which 
can be simply solved by the inversion of a 4 x 4 matrix at 
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each point and for each neutrino species. To ensure that this 
procedure is reasonable, however, we need an estimate of the 
time stepping error - which will help us choose the time step 
dt. We consider four different types of errors. 


when E < 0 and 

7 VFiFj-E 2 

Ej? — ~ ~ 

oiReiE + a A bs max (E) 


(D9) 


• Large fluxes may cause the explicit portion of the time 
stepping algorithm to be unstable. We thus require that 


when FiFj > E 2 . In practice, we find that these 
errors are generally smaller than ee^f- 


apF + o^Abs max ( E ) 
\dj(aFi -poE)\ 


(D4) 


at any point where the denominator is positive. For 
post-merger evolutions, we typically choose op = 0.3 
and c^Abs = 0.001. The maximum is taken over the 
entire computational domain. 


largest error among all points, all neutrino spec ies, and all five 
errors e e , £e , £y , , e' f . If the inequality (D4) is not satisfied 
or if e > 10 , we reject the last time step and start again with a 
new time step given either by ( |D4| ) or by 


dt' = dt 



(D10) 


• Strong coupling of the neutrino radiation to the fluid 
can cause large oscillations in the fluid quantities when 
the radiation transport code is first turned on, causing 
the evolution to be unstable. These swings, when they 
occur, are particularly noticeable in the evolution of the 
electron fraction Y e . Accordingly, we define 


jAnj po 

<^Rel Po + C^Abs max (po) 


(D5) 


where A Y e is the change in Y e over the last neutrino 
time step due to the coupling between the fluid and the 
neutrino radiation. For post-merger evolutions, we use 

<^Rel = 0.01. 


• The implicit portion of the time step might be inaccu¬ 
rate, for example because of the linearization of S a . 
Accordingly, we solve the implicit problem twice: once 
with a time step dt, and once using two time steps 
dt/ 2. The explicit terms are kept identical for both time 
steps, but the linearization of S a is recomputed at the 
intermediate point t + dt/ 2. For any variable U, we 
thus have two estimates U\ and U 2 for U(t + dt), us¬ 
ing respectively one and two intermediate time steps. 
We can then obtain a second-order accurate estimate 
U(t + dt) = 2 U 2 — U\ and an estimate of the error, 
SU = U 2 — U\. We then define 


SE 

a Re \E + a A bs max (E) 
-fVSFiSFj 

a Re \E + a Abs max (E) 


(D 6 ) 

(D7) 


Otherwise, we accept the time step and set the next time step 
for the evolution of neutrinos to be 


dt' = dt 


/ 


0.9 

max (e, 0.3) 


(Dll) 


Finally, if dt' is larger than the time step required for the neu¬ 
trino evolution to catch up with the fluid evolution, we take a 
time step bringing the two sets of equations to the same time. 
In practice, for most of the post-merger evolution, the code 
ends up taking a single neutrino time step for each GR-Hydro 
time step, as can be expected for a quasi-equilibrium radia¬ 
tion field. The adaptive time stepping algorithm is thus mostly 
useful to get stable evolutions when the neutrino radiation is 
first turned on, as well as during the plunge of the neutron 
star into the black hole. Occasionally, taking 2-3 neutrino 
time steps per hydrodynamical time step can also be neces¬ 
sary when time stepping is limited by the Courant condition, 
which is more restrictive for neutrinos than for the fluid. 


Appendix E: Code Tests 

In order to test the ability of our code to perform neu¬ 
trino transport simulations in compact binary mergers, we per¬ 
form a series of tests showing that the moment formalism has 
been properly implemented in SpEC, that it suffers from the 
same known limitations as in other codes, and finally that 
the choices made when averaging emissivities and opacities 
still allow us to reproduce reasonably well the output of an 
energy-dependent code in a spherically symmetric test evolv¬ 
ing a post-bounce configuration taken from a core collapse 
supernova simulation. 


as the errors from the implicit time step. 

• Errors in the evolution of E and Fi can cause E 
to become negative, or the fluxes to violate causality 
( 7 lj FiFj > E 2 ). To be safe, we additionally define 


\E\ 


e E — 


a Re \E + a Abs max (E) 


(D 8 ) 


1. Single beam in flat spacetime 

The simplest test that we perform is the propagation of a 
beam of radiation in vacuum and for a flat spacetime. We use 
a 3-dimensional grid with spacing Ax = 0.01 and bounds 
0<x<l, 0<y< 0.2 and 0 < z < 0.2. The field 
variables are frozen for x < 0.1, with the condition E — 1, 
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FIG. 16. Energy density for the single beam test, after At — 0.4. 
The black arrows show the flux Fi . 


FIG. 18. Energy density for the shadow test after At = 1. The 
arrows show the flux Fi . 



The Ml closure is known to perform well for such a test, 
in which we expect to see the shadow of the optically thick 
sphere in the fields for x > 0.5 (see e.g. ED). Indeed, we 
observe a clear shadow behind the sphere, with the formation 
of two independent beams each having properties similar to 
those observed in the previous tests. Our implementation of 
the Ml formalism thus appears to properly project the shadow 
of optically thick objects. 


FIG. 17. Energy density for the oblique beam test, after At — 0.2. 
The arrows show the flux Fi. 


Fi = (0.9999999,0,0) inside of a “beam” confined to 0.05 < 
y < 0.15, 0.05 < z < 0.15, and are suppressed by a factor of 
10“ 10 in the region outside of the beam. 

The results of this test are shown in Fig. |T6} at a time 
t = to + 0.4. As expected, the beam propagates in a straight 
line and at the speed of light. Numerical errors cause a slight 
widening of the beam, but by only one grid spacing at an 
energy density around 10 _3 F b eam. The beam front is also 
smoothed upstream of the beam, with an exponential decay 
of about 10“ x / 2 per grid point. This pattern establishes itself 
quickly as the beam leaves the frozen region, and then propa¬ 
gates at the speed of light with the beam. 

We then performed the same test, but for a beam no longer 
propagating along a coordinate direction. In this case, we 
freeze the region with x < 0.03 or y > 0.17, and set 
Fi = (0.8, —0.6, 0) in the beam region. The results are shown 
in Fig. [IT] The main difference with the aligned beam is that 
the edges of the beam are not as sharp (while the front of the 
beam is slightly sharper). The exponential decay of the energy 
away from the beam is now about the same in all directions. 


2. Shadow 


We now consider a uniform radiation field with E = 1 
and Fi = (0.9999999,0,0), hitting a sphere of optically 
thick material. In this test, we set the absorption opacity 
K a = 10 6 within a sphere of radius rs = 0.05 and centered 
on cs = (0.5, 0.1, 0.1). The grid is identical to the one used 
for the beam tests. The results of the evolution are shown in 


Fig. 18 


3. Single beam in curved spacetime 


As our simulations require the evolution of the neutrinos 
close to a black hole, we want to determine how well the 
Ml formalism performs in a strongly curved spacetime. To 
do so, we perform another set of tests on beams propagat¬ 
ing in vacuum, but now in a black hole spacetime. These 
tests are largely inspired from those presented in McKinney 
et al. l8P (on their Fig. 13 and Sec. 5.9). The tests in 1(811 
are performed in 2D in spherical polar coordinates, however, 
while the SpEC code always uses 3D cartesian grids for radia¬ 
tion hydrodynamics. For the sake of comparison, we thus first 
consider a “Kerr-like” spacetime in which the metric is set by 

n Kerr ^ K err met _ 


X , y, z) = g* h f rr (0, x, 0, z), with g^ v 


_ n Kerr( 

" nv V 

ric for a non-spinning black hole, in Kerr-Schild coordinates. 
This metric is unphysical, but allows us to consider effectively 
2D beams without having to develop an axisymmetric version 
of our code. In all tests, the beams are emitted from a region 
of width Az = M in the x = 0 plane, in which E = 1 and 
the flux is chosen so that FiF 1 = 0.998F 2 and aF l — (3 l E 
is along the x-axis. We use a grid spacing Ax = 0.2M (i.e. 
the beam is initially 5 grid points wide), where M is the black 
hole mass. 

We first consider a beam emitted from a relatively large dis¬ 


tance, 7 M < z < 8M, shown in Fig. [T9] The energy density 
in the beam decreases as the neutrinos get farther from the 
black hole, in part due to the gravitational redshift (about a 
10% effect over the length of the beam on the grid), and in 
part due to the spreading of the beam caused by the gravita¬ 
tional bending of the radiation. Within the accuracy shown 
in the previous tests, the beam remains otherwise well con¬ 
tained within the region delimited by the null geodesics de¬ 
fined by x(t = 0) = (0,0,8M), dr/dt{t = 0) = 0 and 
x(t = 0) = (0,0,7 M), dr/dt(t = 0) = 0 (shown as solid 
black lines on the figure). 

We now move on to a beam initially closer to the black hole, 
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6dT“ 8.0 


FIG. 19. Energy density in steady-state for a 2D beam in a Kerr- 
like spacetime, with the beam emitted from r = 7M — 8 M and 
M the mass of the non-spinning central black hole. The solid lines 
show null geodesics bounding the expected location of the beam. 
The simulation is performed in full 3D, but the metric and beam are 
independent of z. 


0.5 



4.0 


FIG. 21. Same as Fig. |20] but using the full 3-dimensional Kerr- 
metric and a beam of finite transverse size. 


2.0 


0.0 


FIG. 22. Contours of constant energy density for the crossing beams 
test. The contours are for 10%, 1% and 0.1% of the original energy 
density of the beams. 






FIG. 20. Same as in Fig. [T9] but with the beam emitted from r = 
3 M — 4 M. As before, the solid lines trace null-geodesics. 


with 3 M < z < AM. The inner edge of the beam then lies 
on the unstable circular photon orbit. Both the gravitational 
redshift effect and the divergence of the null geodesics are 
much stronger in this case, so that the beam widens and its 
energy density decreases quickly as the beam orbits the black 
hole. The results of the simulation are shown in Fig. [20| As 
before, we can check that the beam remains mostly within the 
region predicted for free-streaming massless particles. 

Finally, we consider the same configuration but in full 3D. 
The background metric is now that of a non-spinning black 
hole in Kerr-Schild coordinates. In 3D, the null geodesics fol¬ 
lowed by our beam converge towards the equatorial plane, and 
should all cross on the 2 = 0, y = 0 axis. The numerical re¬ 
sults are shown in Fig. 21 As in the 2D test, the widening of 
the beam is consistent with the predictions obtained by trac¬ 
ing null geodesics. The decrease in energy density due to the 
widening of the beam in the xy -plane and to the gravitational 


redshift is however compensated by the collapse of the beam 
along the z-axis, at least until the beam crosses the x-axis. In 
the z < 0 plane, the beam should rapidly widen again along 
the y- axis. In practice, however, the beam remains slightly 
more collimated than expected in that region. Although the 
qualitative differences are fairly mild in this case, it is a first 
indication of the problems intrinsic to the use of the Ml for¬ 
malism, which we will made clearer in the next test. 


4. Crossing beams 

The problems with converging beams, which we alluded to 
in the previous section, can be more easily observed if we con¬ 
sider a much simpler setup: two crossing beams in flat space- 
time. As we do not include any interaction between neutrinos, 
except indirectly through the effective opacity due to the in¬ 
verse reactions of thermal processes, the two beams should 
simply pass through each other. However, this does not occur 
when using the moment formalism with the Ml closure. We 
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FIG. 23. Energy density of neutrinos at t = 2 for the emitting sphere 
tests with k = 5 and k = 10 10 . The numerical results are shown by 
circles/squares while the solid lines are the analytical predictions. 


show in Fig. 22 what happens when such a system is evolved: 
the beams collide and form a wider, single beam along the 
average direction of propagation of the crossing beams. This 
is, effectively, due to the difference between the approximate 
form of the second moment of the distribution function in the 
analytical Ml closure, and its true second moment. 

This inability of the Ml formalism to deal with crossing 
beams, and the inaccuracies existing in the presence of con¬ 
verging beams, means that the evolution of the moments of 
the neutrino distribution function in the region along the polar 
axis of the black hole, as well as close to the inner edge of the 
disk, is not entirely reliable. The Ml formalism performs very 
well for diverging free-streaming neutrinos (which is what we 
see in most of the regions surrounding the disk), and in opti¬ 
cally thick regions (as in the core of the disk, see next test). 
But its performance in regions in which higher moments of 
the distribution function are required to properly model the 
neutrino pressure is, by definition, quite poor (see also ED). 


5. Optically thick radiative sphere 

The tests presented above are mainly intended to determine 
how well our code can evolve the neutrino moments for var¬ 
ious geometric configurations in the free-streaming limit. In 
order to evolve neutron stars and the optically thick accretion 
disks created as a result of neutron star mergers, we also need 
to determine whether we can properly model optically thick 
regions and, more importantly maybe, how closely we can re¬ 
produce the transition between optically thick and optically 
thin regions. Indeed, these will determine how well we can 
predict the neutrino luminosity from neutron stars and accre¬ 
tion disks. 

We first consider optically thick radiative spheres, for 
which we have an analytical solution for the distribution func¬ 
tion (82) : 


f(r, M) = b[l - exp(-Ks(r, g,))} (El) 


for a sphere of radius r s and absorption opacity k,, and for 
fj, = cos($) with 9 the angle between the neutrino momentum 
and the radial direction. The function s(r, /i) is given by 

s{r , g)mrg + r s g(r, g) [r < r s ; -1 < g < 1} , (E2) 

= 2 r s g(r, g) [r > r s ; \j\- (^) < g < 1] , 

= 0 otherwise , 

with 



We can then obtain the exact solution by numerical integra¬ 
tion of /(r, /i) over fi. In this test, we consider two differ¬ 
ent regimes. First we use a sphere of moderate optical depth 
(r s = 1, ft = 5), which is fairly typical of the conditions 
in post-merger accretion disks. Then, we use a sphere of ex¬ 
tremely high optical depth (r s = 1, k = 10 10 ), which pro¬ 
vides a sharp neutrinosphere similar to what may be found 
at the surface of a neutron star. In both cases, we use a 3- 
dimensional Cartesian grid with grid spacing Ax = 0.04. The 
results of the numerical evolutions are compared with the ex¬ 
act solution in Fig. [23] We find good agreement between the 
two solutions. For the most optically thick case, this is mostly 
a consequence of the corrections to high optical depth regions 
described in Appendix |C| The configuration with tz = 5, on 
the other hand, has a very small correction to the fluxes, as 
tzAx = 0.2 < 1. 


6. Neutrino emission and coupling with matter in a 
post-bounce configuration 


As a last test of our code, we consider a ID profile con¬ 
structed as a spherical average from a 2D, Newtonian, multi¬ 
neutrino energy, multi-neutrino angle neutrino transport sim¬ 
ulation at 160 ms after core bounce (83). Since we spheri¬ 
cally average the 2D simulation, do not use precisely the same 
nuclear equation of state and neutrino opacities, and freeze 
the hydrodynamics, the initial profile is slightly out of equi¬ 
librium. Therefore, this provides a good test of the lepton 
and energy coupling and a venue for comparing with energy- 
dependent transport. Similar tests were carried on with this 
profile for testing Monte Carlo transport [[84]. We compare 
the SpEC results with the output of the GRID code (39) . For 
this test, GRID uses an energy-dependent Ml scheme with 
12 energy groups, and was itself shown to agree well with 
the results of full transport codes (85) . During this test, the 
matter density is fixed and the fluid is assumed to be at rest. 
But the internal energy and electron fraction of the fluid are 
coupled to the neutrino evolution. We perform this test with 
various choices of gray schemes: with the standard SpEC 
methods described in Appendix [B] (SpEC-Std), with the av¬ 
erage energy of neutrinos always obtained by assuming equi¬ 
librium between the neutrinos and the fluid (i.e. ignoring the 


correction given by Eq. B8 SpEC -Tfi), and without the cor¬ 
rection [B6] to the charged-current emissivities in low-opacity 
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FIG. 24. Temperature after 8 ms of evolution for the evolution of 
a post-bounce configuration, in the region in which the temperature 
evolution is most affected by the gray approximation. 



FIG. 25. Electron fraction Y e after 8 ms of evolution for the evolution 
of a post-bounce configuration. 


regions (SpEC-? 7 x). The SpEC simulations, which are per¬ 
formed with the full 3D code assuming octant symmetry, have 
a fairly coarse resolution of 6 km at the lowest resolution, 
3 km at the medium resolution, and 1.5 km at the highest res¬ 
olution. 

We first compare radial profiles of the fluid variables T and 
Y e after 8 ms of evolution, shown in Figs. 24||25 The tem¬ 
perature profiles agree well with the GRID results, as long 
as we correct the average energy of neutrinos according to 
Eq. |B8| If we do not, the absorption of neutrinos in the low- 
density regions is widely underestimated - and in particular, 
the code completely misses the existence of a gain region at 
r > 100 km. With correction B8 the heating in the gain re¬ 
gion is reproduced better, but now occurs faster than expected. 
This is due to an overestimate of the average energy of neutri¬ 
nos. Whether or not we modify the emissivities according to 
equation |B6| does not appear to affect the temperature evolu¬ 
tion. 


The evolution of the electron fraction Y e is more sensitive to 
the choices made when energy-integrating. With our standard 
choices, we find results close to GRID, while the composi¬ 
tion is widely inaccurate when neglecting correction |B 8 1 and 
the equilibrium composition can in some regions be wrong by 
A Y e rsj 0.05 when neglecting correction B6 The grid resolu¬ 
tion does not significantly affect these results, except close to 
r = 0. 


Finally, we observe over 20 ms the neutrino luminosity ex¬ 
tracted on a sphere of radius r = 250 km, and compare the 
average energy-squared (e 2 eak ) predicted by the leakage (and 
used to compute opacities in colder regions) and the same 
quantity measured by the energy-dependent Ml code. We find 
that the luminosities are accurate to ~ 10% for v e and v x , and 
to 20% — 30% for v e . For comparison, the luminosities pre¬ 
dicted by the leakage scheme are off by factors of 2-3 for v e 
and v e (but are very good for v x ), while the errors due to the 
finite grid resolution are ~ 5% —10%. The energies yj (e 2 ) are 
typically overestimated by the leakage scheme, by up to 50%. 
This causes the larger-than-expected heating rate observed in 
the gain region (r > 90 km). The reason for these large er¬ 
rors is that the neutrino energies in the hot, highly degenerate 
matter present at high optical depth are very large and sig¬ 
nificantly affect the average energy (due to weighting of the 
electron number emission rate by the square of the neutrino 
energies). In practice, however, these high energy neutrinos 
are thermalized as they propagate through the optically thick 
regions. Significant improvements in y/ (e 2 ) can be obtained 
if we define 


<e 2 > = 


/-Mx)(e^(x)) min (1, exp ^threshold -T])dV 
f R v (x) min (1, exp[r t hreshoid - t}) dV 


(E4) 

where R v is the neutrino number emission predicted by the 
leakage scheme, and r is the optical depth. Any value of the 
threshold Threshold ^ 1 — 10 allows us to recover yj (e 2 ) to 
within 20% for v e and z7 e , and within 30% for v x . In all cases, 
the finite resolution error is negligible compared with the er¬ 
ror due to the gray approximation. We have checked that these 
variations are indeed due to the gray approximation, and not 
to the specific implementation of the Ml scheme in SpEC: an 
energy-dependent version of the Ml code in SpEC, which can 
easily be used in this case due to the relatively low compu¬ 
tational cost of the test problem and the lack of gravitational 
redshift or velocity gradients, agrees with the GRID results. 

These results are already much better than if we got the 
neutrino energies by assuming equilibrium with the fluid: in 
the gain region, the equilibrium average neutrino energies are 
y/ (e 2 ) ^5 — 10 MeV while the actual average neutrino en¬ 
ergies are yj (e 2 ) ~ 10 — 25 MeV. And indeed, with this 
corrected (e 2 ), heating in the gain region agrees well with the 
results of the GRID simulations. The correction also causes 
changes to the evolution of Y e , with results with and without 
correction |B6| to the emissivities now bracketing the GRID 
results (with typical errors A Y e ~ 0.02). The simulation 
without correction B6 still slightly overestimates Y e around 
r ^ 100 km, while with correction |B6| the simulation now 
slightly underestimates Y e in that same region. 
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Such errors would be very significant in the context of 
core collapse supernova simulations, where properly model¬ 
ing neutrinos is critical to the explosion mechanism. In com¬ 
pact binary mergers, however, the impact of neutrinos is more 
modest. The computation of (e 2 ) from the leakage scheme is 
also more accurate in post-merger accretion disks: the tem¬ 
perature does not vary as much within post-merger accretion 
disks as in the test presented here, the disks are only mod¬ 
erately optically thick (r < 10), and matter within the disks 
reaches densities of at most po ^ 10 12 g/cm 3 , at which only 
electrons are degenerate. Using emissivities with and without 
correction |B6[ we can get an estimate of the errors due to the 
use of a gray scheme. These errors, discussed in the main text 
of the paper, are much smaller than the errors observed in this 
test. Using the Ml formalism instead of our leakage scheme 
then significantly improves our ability to determine the cool¬ 
ing time and composition evolution of the disk, and gives us a 
reasonable first estimate of the energy deposition in the corona 
(except along the polar axis of the black hole, where the Ml 
approximation is unreliable). 


Appendix F: Accuracy of the post-merger evolution 


Most of the discussion of our simulations of the formation 
of an accretion disk after a neutron star-black hole merger 


(Sec. IV) focuses on the qualitative features of the system, and 
on the differences between various algorithms for the evolu¬ 
tion of neutrinos. In this section, we will discuss the accu¬ 
racy of these results, and argue that the features of the system 
discussed in Sec. [TV] are appropriately resolved by our simu¬ 
lations. We identify four main sources of error. The first is 
due to the finite resolution of our numerical grid during the 
simulations. The second is the numerical error in the simula¬ 
tion of the neutron star-black hole merger before we turn on 
the neutrino transport code, which causes errors in the initial 
conditions used for our simulations. The third is the approx¬ 
imate treatment of the neutrinos, and in particular the use of 
a gray scheme and of the Ml closure. And the fourth comes 
from only turning on the neutrino transport and changing the 
treatment of low-density material about 6.1 ms after merger. 
The first two can be estimated through the use of simulations 
at lower and/or higher resolution. The last two are more dif¬ 
ficult to assess, and will only be rigorously measured through 
the use of more advanced (and more costly) simulations. Un¬ 
til such simulations are available, we have to rely on simpler 
estimates of these errors. 

To determine the importance of the first source of error, we 
performed simulations of a post-merger accretion disk using 
60 3 and 140 3 points for each level of refinement, instead of 
the 100 3 points used in our ’’standard” configuration. The 
higher resolution simulation was only evolved for 1.1ms, to 
verify that the solution was converging with resolution. Even 
the simulation using the coarsest grid shows surprisingly good 
agreement with our standard runs. The average temperature in 
the disk agrees to A (T) ^0.2 MeV and the electron fraction 
to A(y e ) ~ 0.01. The neutrino luminosities agree to better 
than 20%, and the average neutrino energies within 0.5 MeV. 


This is comparable to the differences observed between var¬ 
ious choices of gray approximations, and much smaller than 
the differences between the leakage and Ml simulations. Ad¬ 
ditionally, the difference between the standard and high reso¬ 
lutions is, for all observed quantities, more than a factor of two 
smaller than the difference between the low and standard res¬ 
olutions. Most of those differences arise immediately after the 
neutrino transport is turned on, while evolution at later times is 
very similar for all numerical resolutions. Accordingly, we do 
not expect numerical resolution during the post-merger evolu¬ 
tion to be a significant source of error. 

The numerical error in the simulations before we turn on 
neutrino transport can easily be determined from the lower 
resolution simulations of the same system performed in Fou- 
cart et al. 2014 l24l . The largest error in these simulations 
was in the determination of the properties of the dynamical 
ejecta, which does not concern us here as that material is al¬ 
lowed to escape the numerical grid. At the time at which we 
begin the simulations with neutrino transport, we otherwise 
find relative errors of less than 10% in the total mass outside 
of the black hole, average temperature and average electron 
fraction in the disk. Considering that the grid spacing used 
in lf24l was about a factor of two coarser than the grid spacing 
used in the standard simulations of this paper, these are clearly 
overestimates of the numerical error. The only numerical error 
which could affect our results is thus the initial temperature 
of the disk ,which generally decreases with resolution. The 
disk could be, on average, about 0.4 MeV cooler at infinite 
resolution. This is only slightly smaller than the difference 
between leakage and Ml simulations, but does not affect the 
conclusion of this paper: both the leakage and Ml simulations 
presented here use the same imperfect initial data. We should 
also note that even our imperfect initial data remains a much 
better starting point than the commonly used alternative, an 
equilibrium torus of constant entropy. 

The effect of an approximate treatment of the neutrinos is 
discussed in Sec. |IV G[ to which we refer the interested reader. 
Here we will simply note that our rough estimate of the er¬ 
ror due to the use of a gray scheme is comparable to our 
estimate of the error due to numerical resolution, and much 
smaller than the differences between simulations using the 
Ml scheme and the leakage scheme. However, a true mea¬ 
sure of the error would require a simulation using an energy- 
dependent transport scheme without the weaknesses of the Ml 
closure. Such a simulation is currently out of reach for our 
code. 

Finally, we have to consider the impact of suddenly turn¬ 
ing on the neutrino transport code 6.1 ms after merger, and of 
the use of a higher threshold for the application of atmosphere 
corrections at earlier times. The only way to rigorously mea¬ 
sure this would be to perform a simulation using the Ml code 
and a lower atmosphere threshold, starting before the disrup¬ 
tion of the neutron star. Now that we have a well-tested Ml 
code and good estimates of the atmosphere thresholds which 
should be used after merger, we intend to perform such sim¬ 
ulations. However, the disruption of the neutron star and ac¬ 
cretion of the neutron star core onto the black hole is the most 
computationally intensive part of the evolution, while we be- 
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lieve that the simulations performed here capture the most 
important properties of the merger and post-merger evolution 
(except for the absence of magnetic fields). Indeed, we begin 
the simulation around the time at which the neutrino luminos¬ 
ity and the electron fraction begin to significantly increase due 
to the heating of the disk and the emission of neutrinos. We 
also find that the disk outflows are not immediately launched 
when we decrease the atmosphere threshold. At the beginning 
of the simulation, those outflows are still contained by mate¬ 
rial from the tidal tail falling onto the forming disk. It is only 
about 2 ms later that outflows begin to form above the disk. A 
potentially more important source of error is the transient oc¬ 


curring when the transport code is turned on. To study its ef¬ 
fects, we consider two possible initializations of the moments 
of the neutrino distribution function. First, we initialize them 
assuming that the neutrinos are in equilibrium with the fluid, 
in which case the neutrino energy density is overestimated, 
and the initial transient consists in a transfer of energy from 
the neutrinos to the fluid and a slight decrease of the average 
electron fraction. Then, we initialize the neutrino energy den¬ 
sity to a negligible value, in which case the transient has the 
opposite effect. We verify that after about 2 ms the two solu¬ 
tions are in good agreement (at early times, this transient is of 
course the largest source of error in the simulation). 



